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Abstract 


The  modeling  and  simulation  work  in  this  grant  was  directed  by  Prof.  D.  Scott 
Stewart,  (the  Stewart-Group)  and  the  experimental  work  and  data  collection  was 
directed  by  Prof.  Nick  Glumac  (the  Glumac-group).  Professors  Stewart  and 
Glumac  and  Drs.  John  Bdzil  and  Joseph  C.  Foster,  collaborated  on  experimental 
designs,  that  benefited  from  the  concurrent  theoretical  and  computational 
modeling.  The  effort  described  in  this  report  fully  integrated  modeling  and 
experiments  to  study  the  energy  release  processes  that  thermitic  and/or 
exothermic  intermetallic  reactive  materials  experience  when  they  are  subjected 
to  sustained  shock  loading. 
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I.  Objectives 

An  objective  of  this  grant  was  to  integrate  modeling  and  experiments  to  develop 
an  increased  scientific  understanding  of  the  condensed  phase  energy  release 
processes  in  the  microstructure,  that  thermitic  and/or  exothermic  intermetallic 
materials  experience  when  they  are  subjected  to  sustained  shock  loading.  Data 
from  highly  spatially  and  temporally  resolved  experiments,  supported  by  our 
modeling  and  simulation  efforts  were  directed  at  building  and  validating  a  model 
that  describes  the  basic  physical  processes,  including  solid  state,  multi-phase 
transport  and  reaction/diffusion  mechanisms.  Such  experimentally  validated 
models  are  invaluable  for  the  design  of  the  intended  applications.  Our  efforts 
provided  a  unique,  multi-disciplinary  research  experience  for  graduate  students 
and  young  scientists  affiliated  with  and  supported  by  this  effort. 

I.  1  Background 

The  Department  of  Defense  (DoD)  and  the  Defense  Threat  Reduction  Agency 
(DTRA)  has  had  a  long-  standing  interests  in  designing  reactive  material  (RM) 
whose  energy  release  is  triggered  by  shock/shear  stimulus  and  which  contributes 
to  long  duration,  high-thermal  pulse.  A  main  application  of  thermitic  reactive 
materials  is  the  defeat  of  biological  agents  that  must  be  destroyed  with  an 
extremely  high  temperature  deflagration,  that  is  highly  contained  in  localized  area. 
There  is  a  need  to  avoid  wide-scale  dispersal  of  the  biological  threat  agent  and 
thermal  destruction  is  a  preferred  mechanism.  Materials  that  lie  in  a  group  of 
intermetallic,  and  thermitic  mixtures  are  promising  because  when  the  unreacted 
(solid)  materials  are  combined  to  make  products  materials,  they  often  have  very 
high  (enthalpic)  energy  release  and  often  make  make  mainly  solid  products  so 
that  gas  generation  is  somewhat  minimized.  Our  original  proposal  identified  a 
number  of  target  systems  of  interest  that  included,  the  Titanium-Boron  system 
(metal/intermetallic),  the  aluminum  copper  oxide,  (metal/metal  oxide)  system  as 
well  as  a  number  of  similar  combinations  with  metal/metal  oxide/intermetallic 
components. 

II.  Approach 

In  design  applications,  the  RM’s  are  typically  manufactured  composites,  layered 
or  particulate,  and  the  reactants  are  not  molecularly  pre-mixed.  Therefore  a 
distinct  microstructure  comprised  of  regions  of  separated  reactants,  is  inherited 
from  the  constituents  and  fabrication.  Reaction  of  the  initial  constituents  must 
take  place  at  interfaces  between  molecularly  pure  materials.  The  initial  shock 
passage  serves  as  a  trigger  that  concentrates  energy  by  collapse/compaction 
and  yield  of  nominally  solid  materials,  that  surely  experience  temperature  rise 
through  dissipative  yield  mechanisms.  But  immediately  and  subsequently  the 
overall  stress  level  drops  to  near  ambient  conditions  of  the  larger  surroundings, 
while  high  temperature  remain  in  localized  regions  at  the  interfaces  between 
materials,  The  post-shock  energy  release  that  is  related  to  chemical  reaction  is 
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controlled  entirely  by  anerobic,  condensed  phase  processes  in  the  RM.  The 
project  scope  was  based  on  a  scenario  from  application,  where  a  shock  is  driven 
into  RM  by  a  high  explosive,  and  the  subsequent  reaction/combustion  of  the 
materials  is  observed  and  modeled.  Below  Figure  6.  from  the  original  proposal 
provides  an  excellent  conceptual  schematic  of  this  effort 


Figure  6.  A  schematic  depiction  of  the  normal  shock  loading  and  shear  loading  experiments 
described  in  Subsection  A)  and  B),  respectively.  The  bottom  of  the  figure  shows  a  seguence  in  an 
x,t-space,  with  two  reactant  materials  A  and  B  subjected  to  a  driving  shock  that  passes  through 
material  A  and  then  transmits  the  shock  into  B,  with  a  reflected  shock  traveling  back  through 
material  A.  If  the  materials  are  pure  reactant  constituents,  then  reaction  first  starts  at  the  material 
interface.  If  the  reaction  A  +  B  ->  C  proceeds  and  is  sustained,  then  diffusion  of  A  and  B  into  the 
diffuse  region  of  product  C  maintains  the  reaction  near  the  interface.  This  simple  diffusion-reaction 
model  is  a  stand  in  for  the  more  realistic  thermite  and  intermetallic  reactions. 


Figure.  1.  A  reproduction  of  Fig.  6  from  the  original  proposal  that  describes  the  basic 
scenario  and  premise. 

At  the  time  of  the  inception  of  this  proposal  (2009)  even  to  the  time  of  this 
summary  report  (2015),  with  a  few  exceptions  that  include  work  done  by  our 
combined  groups,  there  have  been  few  high-resolution  experiments  that  attempt 
to  resolve  the  energy  release  processes  in  the  microstructure  and  few  models  of 
the  condensed  phase  processes  that  are  truly  predictive.  While  experiments  are 
limited  by  optical  issues,  modeling  and  simulation  issues  have  been  limited  by 
the  necessary  requirements  to  simultaneously  model  deformation,  reaction  and 
diffusion  in  condensed  materials.  Our  efforts  focused  on  exactly  these  issues. 

Two  major  sets  of  accomplishments  were  made  during  our  efforts.  1)  On  the 
experimental  side,  the  Glumac  group  executed  a  set  of  experiments,  [25]  on 
shock  compaction  of  RM  that  was  able  to  visualize  the  characters  of  ignition  and 
propagation  of  reaction  in  thermite  and  intermetallic  systems.  2)  On  the 
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experimental  side,  the  Stewart-group  was  able  to  develop  new  models,  (e.g.  the 
Gibbs  formulation  [18])  that  is  fundamentally  based  on  the  concept  of  stress  and 
temperature  equilibrium  at  the  smallest  continuum  scales  at  all  points  in  materials 
mixtures.  Hence  stress  and  temperature  are  fundamental  state  variables,  that 
define  Gibbs  potentials  for  the  materials.  This  new  approach  allowed  our  group  to 
develop  new,  material  models  for  reaction  and  diffusion  and  deformation  in 
condensed  phase  reactants.  Both  of  these  main  accomplishments  have  been  (or 
are  currently  being)  documented  in  archival  publications  by  the  members  of  our 
combined  group,  [1 5, 1 6, 1 7, 1 8, 1 9,20,21 ,24]  These  upcoming  and  planned 
archival  publications  will  cite  DTRA  support,  [25,26,27,28] 

What  follows  below  is  a  brief  overview  description  of  the  "Striker  experiment" 
carried  out  in  the  Glumac  lab,  and  a  summary  collection  of  the  experimental 
results  obtained  so  far.  This  is  followed  by  a  brief  description  of  some  of  the  key 
modeling  and  simulation  results.  These  are  the  main  achievements  of  this  grant. 
A  complete  description  of  this  work  and  its  chronological  development  is  found  in 
Section  VII,  which  is  a  concatenation  of  our  annual  reports. 

III.  Synopis  of  Work  Accomplished 

III.  1  Experimental  Work 

The  first  experimental  task  was  to  define  an  representative  experiment,  that 
would  capture  key  features  of  the  ignition  of  a  reactive  material  (RM)  that  was 
subjected  to  shock  and  rapid  compression.  Our  entire  team  of  senior 
investigators,  experimentalists  Glumac  and  Foster,  and  modelers  Stewart  and 
Bdzil,  worked  together  to  come  up  with  a  conceptual  design  we  called  the  "Striker 
experiment".  In  this  experiment  a  fabricated  RM  sample  that  was  approximately  a 
3/8  inch  cube,  was  inserted  into  a  steel  confinement  cell.  On  either  side,  steel 
slider  "striker  bars"  were  inserted  and  tamped  against  the  RM  sample.  Then  on 
each  striker  bar,  detonators  (such  as  Reynolds  RP-80  detonators)  were  placed 
on  the  striker  bar  in  order  to  generate  a  shock  wave  and  push  the  striker  bars 
into  the  sample.  A  PMMA  window  for  optical  viewing  was  placed  on  one  side  of 
the  confinement  cell.  On  firing  of  the  detonators,  reactive  events  at  the  boundary 
of  the  RM  and  the  viewing  window  were  imaged  in  a  1  cm  square  view  area  by 
means  of  high  magnification  lenses,  with  very  high  speed  cameras.  Different 
classes  of  RM  materials  were  subjected  to  this  experiment.  Figures  1  and  2  (this 
report)  shows  an  overview  of  the  experiment,  the  basic  setup  (circa  2012),  and 
the  imaging  setup  up  for  the  second  generation  experiments. 
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Overview  of  the  Experiment 


•  Explosively  driven 
pushers  into  RM 
sample.  Watch  with 
HSC, 


•  Can  achieve  grain 
scale  spatial  resolution 
and  sub-us  time 
resolution 


Figure  1.  Overview  of  the  Striker  experiment.  Striker  bars  driven  by  detonators  drive  compression 
bars  into  the  RM  sample  and  direct  imaging  of  burning  events  is  made  through  an  optical  window. 


Basic  Setup 


Figure  2.  Basic  set  up  of  the  Striker  experiment.  The  detonators  are  embedded  in  fragments 
shields  to  prevent  excessive  damage  to  the  confinement  cell,  that  was  used  in  multiple 
experiments. 
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Experimental  Methods: 
Imaging 


Sample  Holder 
and  Base 


Reversed 
Telephoto  Lens 


Phantom  5 
High  Speed  Camera 


Reactive  sample  is  imaged  usini 
80-200  mm  reversed  telephoto 
lens 

Optical  shields  reduce  intense 
emission  from  detonators 

Phantom  camera  records 
images  at  10,000  or  19,000  fps 

Full  sample  and  magnified  field 
of  view  tested  for  each  material 


Calibration 

Grid 


Spatial  calibration  accomplished 
using  calibration  grid 


Figure  3.  Basic  imaging  set  up  of  the  Striker  experiment.  A  reversed  telephoto  lense,  with 
intensity  filters  and  high-speed  phantom  camera  were  key  features  of  the  imaging  system. 


Tests  were  carried  out  on  a  suite  of  RMs  that  were  manufactured  in  the  Glumac- 
lab.  The  basic  materials  were  made  from  mixtures  of  metal  and  metal  oxide 
powders  (such  as  aluminum  and  copper  oxide)  to  make  a  thermite  RM  or  a  metal 
and  intermetallic  materials  (such  as  titanium  and  silicon)  to  make  a 
metal/intermetallic  RM.  These  were  combined  and  pressed  to  about  80  %  of 
their  theoretical  maximum  density.  Titanium/Boron  RMs  were  also  tested  but 
were  found  hard  to  ignite  in  the  Striker  experiement.  However  they  were  fired 
with  mulitple  detonators,  in  project  funded  by  a  grant  funded  supported  by  the  Air 
Force  Research  Laboratory  for  bio-cidal  experiments,  [9,14], 

Figures  4  and  5  show  representative  imaging  results  obtained  from  the  Striker 
experiment,  for  a  stoichiometric  thermite  mixture  of  aluminum  and  copper  oxide. 
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Results:  Stoichiometric  AI/CuO  Thermite 


■  Imaging  and  video  sequences  show  progression  of  the  reaction 

■  Initial  densification  can  be  observed 

■  Propagation  and  flame  speed  can  be  examined 


60.75  ms  81.55  ms  121.35  ms  141.25  ms 


Figure  4.  Images  from  the  AI/CuO  experiment,  with  a  80  %  TMD,  stoichiometric  mixture. 


Figure  5.  Close  up  images,  for  the  same  experiment  shown  in  Fig  4.,  (with  times  not  marked). 
Images  show  evidence  of  a  front  that  follows  compaction  density  contours,  and  secondary  and 
tertiary  reaction. 
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In  addition  to  the  Striker  experiments  shown  above,  another  series  of  experiment 
that  used  a  pure  (fuzed)  quartz,  cylindrical  confinement  cell  was  used.  This 
allowed  greater  optical  access  than  the  steel  confinement  cell.  But  the 
confinement  was  inadequate,  when  the  same  detonator  and  (cylindrical  rod) 
striker  bars  were  used.  The  confinment  came  apart  prematurely.  So  this  set  of 
experiments  while  promising  were  only  partly  successful.  This  redesign  was 
attractive  from  a  modeling  and  simulation  point  of  view.  But  a  revised  version  of 
this  experiment  was  not  attempted. 

I  1 1.2  Modeling  and  Simulation  Work 

Simulations  that  Supported  Experimental  Design 

The  first  efforts  of  our  group  were  to  provide  basic  support  for  the  experiments  in 
the  Glumac  lab.  Simulation  of  the  compaction  events  were  carried  out  with  the 
multimaterial  hydrocode,  ALE3D  (distributed  by  Lawrence  Livermore  National 
Laboratory).  Two  and  three  dimensional  simulations  of  shock/compaction  event 
were  carried  out  in  representative  experimental  geometries  in  order  to  estimate 
timing  and  to  evaluate  the  effects  of  mixing  between  particle  during  the 
compaction  effects. 

A  major  conclusions  was  that  the  effect  of  shock  compaction  is  to  increase  the 
area  where  reaction  can  take  place  at  the  interfaces  between  the  initially 
separated  reactants.  Also  we  found  that  one  must  include  material  strength  in  the 
shock  actuated  mixing,  since  otherwise  excessive  reaction  surface  is  created  in 
the  material  simulations  that  is  not  physical. 

Another  major  conclusion  was  that  the  light/dark  boundaries  observed  in  the 
Striker  experiment  follows  contours  of  densification  (compaction).  The 
densification  process  is  endothermic  where  mechanical  energy  is  concentrated 
into  small  region  of  the  microstructure  and  where  mechanical  dissipation  there 
results  in  large  amount  of  localized  heating.  The  compaction  process  is  realized 
by  the  propagation  of  a  front  and  the  observed  light/dark  front  propagates  in  a 
sequenced  manner  that  follows  the  compaction  front.  Hence  we  found  significant 
evidence  that  burning  takes  place  at  the  separated  reactants  boundaries,  and 
that  the  ignition  of  the  RM  enabled  by  the  compaction  processes. 

Development  of  a  Multi-Component  Reaction  Diffusion  Model  Framework 

One  of  the  main  tasks,  outlined  in  the  original  proposal  was  to  model  multi- 
component  reaction  diffusion  flames  in  condensed  media.  This  led  to  an  effort  to 
identify  thermodynamically  consistent  approaches  and  revisit  the  literature. 
Significant  efforts  were  made  based  on  earlier  work  of  Stewart,  and  later  revised 
to  generate/invent  what  we  now  refer  to  as  the  Gibbs  formulation,  [18].  A  single 
stress  tensor,  and  a  single  temperature  is  assumed  for  the  mixture  with  specified 
Gibbs  potentials  for  all  relevant  species,  and  interaction  energies.  This  work  was 
of  a  fundamental  nature,  within  non-equilibrium  continuum  thermodynamics. 
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The  Gibbs  formulation  model  enjoyed  additional  support  through  grants  from  the 
Office  of  Naval  Research  and  the  Air  Force  Office  of  Scientific  Research. 


Figure  7:  Shows  a  video  frame  from  a  stoichiometric  Ti/Si  experiment 


Figure  8  Simulation  results  that  show  "square"  shapes  apparent  in  the  color 
densification  contours  superimposed  with  a  Schlieren  image,  of  the  Ti/B  mixture 
simulations  the  shape  of  the  initial  compaction  front  contours  and  the  shape  of  the  flame 
front  in  the  experiment. 

Figure  6.  Figures  7  and  8  extracted  from  the  2013  annual  report  that  show  a  qualitative 
comparison  of  the  shape  of  the  initial  compaction  front  contours  from  simulation  and 
and  the  shape  of  the  flame  front  in  the  experiment.  This  comparison  strongly  suggests 
that  in  the  experiments,  the  lead  reaction  front  propagates  along  the  compaction  front 
contours. 
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Modeling  Titanium  Boron  Condensed  Phase  Diffusion  Flames 

Because  of  interest  in  biocides,  the  original  proposal  discussed  in  detail  the 
Titanium/Boron  system.  So  our  initial  examples  of  reactive  flow  modeling 
considered  Titanium/Boron  as  a  representative  (but  significant)  example.  We 
note  that  subsequently  we  found  in  the  experimental  program  that  Ti/B  RM 
mixtures  were  fairly  hard  to  ignite  and  that  systems  like  AI/CuO  or  Ti/Si  were 
more  reactive  and  thus  provided  lots  of  data  in  the  Striker  experiment  setups. 
Nonetheless,  we  note  another  advantage  of  the  choice  of  the  Ti/B  system  for 
fundamental  study,  is  that  the  basic  overall  reaction  is  simple  with  Ti  +  2  B  -> 
TiB2.  Thus  only  three  components  were  needed  for  our  earliest  models  and  first 
treatments  of  condensed  phase  diffusion  flames.  A  version  of  the  Gibbs 
formulation  was  used  that  took  advantage  of  a  simplifying  assumption  that  the 
components  were  fluid-like,  and  that  a  steady  diffusion  flame  that  situated  itself  in 
two  opposed  stream  of  separated  reactants,  was  studied  in  detail  and  recently 
published  [17]. 


Fig.  1.  Schematic  of  the  condensed  phase  counterflow  diffusion  flame,  with  titanium 
entering  from  the  left  and  boron  from  the  right.  The  color  shades  clearly  delineate  the 
planar  mixing  region  and  reaction  zone. 


Figure  7.  Figure  1  extracted  from  a  recently  published  paper  (in  press), "  Diffusion 
flames  in  condensed-phase  energetic  materials:  Application  to  Titanium-Boron 
combustion",  Combustion  and  Flame  (2015) 
http://dx.doi.Org/10.1016/j.combustflame.2015.08.023 
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Modeling  of  the  Ignition  of  an  Aluminum  Copper  Oxide  Thermite 

The  Glumac  experiments  on  aluminum/copper  oxide  thermite  showed  ignition 
and  extinction  behavior.  Evidence  suggested  that  this  takes  place  at  the  interface 
between  the  initially  separated  reactants. 
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FIGURE  1.  The  entire  domain  from  the  viewing  win¬ 
dow  (right-top)  of  the  Glumac  experiment.  The  cutout 
of  image  shows  ignition/extinction  event  in  the  red  cir- 
cle(left),  image  intensity  of  bright  spot(right-bottom) 


Figure  8.  Figure  1  extracted  from  a  recent  APS  proceeding  paper,  that  captures  a  sub¬ 
millisecond  ignition/extintion  event  that  takes  place  at  the  scale  of  the  microstructure  of 
the  aluminum  copper  oxide  thermite. 

Efforts  in  the  Stewart-group  expanded  the  Gibbs  formulation  to  model  the 
features  of  ignition  shown  in  Fig.  8.  This  modeling  effort  considered  of  the  basic 
processes  of  a  thermite,  and  included  no  less  than  eight  components,  and  the 
melting  of  solid  phases  to  liquid.  This  substantially  expanded  the  analysis  of 
multi-component  diffusion.  The  component  considered  included,  minimally  Al- 
solid  and  liquid,  CuO-solid  and  liquid,  Cu-solid  and  liquid  and  Al2  O3  -liquid  and 
solid,  [21,25,27], 
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The  characteristics  of  a  steady  diffusion  flame  that  arises  at  the  interfaces  of  two  condensed  phase  reactant 
streams  that  form  an  opposed  counterflow  are  discussed.  We  assume  that  the  flow  is  due  to  deformation  from 
compaction  or  local  heating  and  thermal  expansion  processes  in  the  microscale  environment  of  composite 
energetic  materials.  As  a  representative  example  of  high  temperature  combustion  of  metal/intermetallic  re¬ 
actants,  the  overall  reaction  of  titanium  and  boron  to  create  titanium  diboride  products  is  considered  under 
near  isobaric  conditions.  The  multi-component  diffusion  description  uses  a  generalized  Fick  formulation  with 
coefficients  related  to  the  binary  diffusivities  defined  in  the  Maxwell-Stefan  relations.  A  fairly  simple  deple¬ 
tion  form  with  Arrhenius  temperature  dependent  coefficients  is  used  to  describe  the  reaction  rate.  Several 
types  of  analyses  are  carried  out  at  increasing  levels  of  complexities:  an  asymptotic  analysis  valid  in  the  limit 
of  low  strain  rates  (high  residence  time  in  the  reaction  zone),  a  constant  mixture  density  assumption  that 
simplifies  the  flow  description,  diffusion  models  with  equal  and  unequal  molecular  weights  for  the  vari¬ 
ous  species,  and  a  full  numerical  study  for  finite  rate  chemistry,  composition-dependent  density  and  strain 
rates  extending  from  low  to  moderate  values.  All  are  found  to  agree  remarkably  well  in  describing  the  flame 
structure,  the  flame  temperature  and  the  degree  of  incomplete  combustion.  Of  particular  importance  is  the 
determination  of  a  critical  strain  rate  beyond  which  steady  burning  may  no  longer  be  observed.  The  analy¬ 
sis  has  a  general  character  and  can  be  applied  to  other  condensed  phase  energetic  material  systems,  where 
reaction  and  diffusion  occur  in  the  presence  of  flow  and  material  deformation. 
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1.  Introduction 

Energetic  materials  are  a  broad  class  of  manufactured  materials 
that  traditionally  comprise  both  propellants  and  explosives,  but  in¬ 
clude  thermite  and  intermetallic/metal  mixtures  as  well.  They  consist 
of  molecular  explosives  and  molecular  oxidizer  crystallites,  like  HMX 
and  Ammonium  perchlorate  (AP),  metal  elements  like  aluminum  (Al) 
and  titanium  (Ti),  metal  oxides  like  iron  oxide,  copper  oxide,  inter- 
metallic  elements  like  carbon,  silicon  and  boron,  and  plastic  binders 
and  resins  like  HTPB.  The  mixtures  are  made  from  elements,  or  a  com¬ 
posite  of  components  that  have  been  processed  and  have  imperfec¬ 
tions  and  contaminants,  such  as  cracks  or  inclusions,  or  have  been 
subjected  to  surface  oxidation,  as  in  propellants.  The  components  in 
the  mixture  prior  to  the  composite  assembly,  all  have  their  own  indi¬ 
vidual  mechanical  and  thermo-chemical  identity. 
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Explosive  compounds  are  pre-mixed  at  the  molecular  level, 
whereas  the  mixtures  of  metals,  metal  oxides  and  intermetallics 
compound  are  initially  separated  constituents  that  react  in  the  vicin¬ 
ity  of  the  initial  material  boundaries,  enabled  by  diffusion  of  reactants 
and  products  through  the  molecular  condensed  phase  regions.  Typi¬ 
cally,  the  individual  components  are  combined  to  make  an  agglomer¬ 
ated  composite  mixture  of  powders  and  other  ingredients.  The  pow¬ 
ders  have  a  particle  size  distribution,  with  dimension  that  varies  from 
one  to  several  hundreds  microns,  and  characteristic  particle  morphol¬ 
ogy.  The  mixtures  are  pressed  or  cast  into  a  mold  for  applications. 

Reactive  energetic  materials,  as  defined  in  a  2004  US  National 
Academy  of  Sciences  report  [1],  generally  combine  two  or  more  in¬ 
ert  solids,  and  an  ignition  source  is  required  to  start  the  chemical 
reaction  at  the  interfaces  of  the  initially  separated  reactants.  Agglom¬ 
erated  composites  can  be  ignited  by  shock  compression  if  the  com¬ 
ponents  comprising  the  initial  reactive  mixture  are  not  highly  dense. 
The  conditions  required  for  ignition  occur  through  the  collapse  of  in¬ 
terstitial  voids,  typically  of  initial  volume  fraction  in  the  range  of  1- 
20%.  During  the  void  collapse,  thermo-mechanical  flows  and  mate¬ 
rial  displacements  occur  due  to  the  densification  caused  by  the  shock 
confinement,  and  lead  to  large  strain  rates  in  the  reactive  material 
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components.  The  resulting  deformations  are  associated  with  energy 
dissipation  that  cause  localized  heating  (hot  spots)  which,  in  turn, 
generates  thermal  events  that  are  needed  for  melting  the  constituent 
components  and  to  initiate  the  reaction  at  the  material  interfaces.  The 
localized  thermal  events  in  the  mixture  also  lead  to  local  streaming 
motions  and  greater  species  mobility  of  reactants. 

Energetic  materials  fall  generally  into  two  broad  classes  based  on 
how  they  are  manufactured:  (i)  agglomerated  composite  mixtures  of 
powders  and  other  reactive  components,  and  (ii)  finely-spaced,  struc¬ 
tured  or  layered  composites  with  specified  interstitial  spacing  be¬ 
tween  the  reactive  components.  The  current  state  of  experimental 
investigations  and  modeling  of  finely-spaced  multilayered  reactive 
materials  has  been  recently  reviewed  by  Rogachev  [2],  Weihs  [3]  and 
Adams  [4],  A  typical  illustration  is  described  by  Weihs  [5],  where  ap¬ 
proximately  10  micron  layers  of  foil  (of  say,  nickel  and  aluminum)  are 
pressed  into  a  laminate  ply.  The  formation  of  nickel/aluminum  prod¬ 
ucts,  for  example,  starts  with  thermal  initiation  (heating)  of  the  com¬ 
posite  material  that  first  melts  the  aluminum  and  allows  the  chemi¬ 
cal  reaction  to  proceed.  The  modeling  effort  focuses  on  the  effect  of 
the  laminate  bilayer  thickness  on  ignition  and  on  the  velocity  of  the 
self-propagating  reactions  of  a  number  of  reactive  material  pairs.  In 
general,  steady  self-propagation  of  a  reacting  front  in  the  laminate 
material  moves  perpendicular  to  the  normal  of  the  planes  of  lami¬ 
nate  plies.  The  speed  of  the  steady  self-propagation  reaction  front  is 
found  to  be  dependent  on  the  ply  spacing  and  composition  [6],  For 
example,  Sraj  [7]  studied  the  response  of  Ni/Al  multilayered  compos¬ 
ites  to  shock  compression;  he  used  the  Sandia  CTH  hydrocode  to  sim¬ 
ulate  the  mechanical  deformation  of  layered  Ni/Al  material  to  shock 
compression  and  investigate  the  effects  of  the  bilayer  thickness  and 
the  shock  velocity  and  orientation  on  the  combustion  process.  When 
the  normal  to  the  plane  of  the  impact  that  creates  the  shock  is  in  the 
same  direction  as  the  normal  to  the  plane  of  the  laminate  plies,  the 
propagation  process  is  found  to  be  highly  unsteady. 

An  alternative  to  the  well-structured  layered  composites  are  ag¬ 
glomerated  materials  made  of  a  mixture  of  reactive  and  inert  com¬ 
ponents.  Glumac  et  al.  [8]  have  recently  reported  results  of  shock 
compaction  experiments  on  porous  materials  that  are  initially  com¬ 
posed  of  two  reactive  components.  Systems  that  have  been  studied 
include  the  aluminum  (Al),  copper  oxide  (CuO)  thermite,  and  the 
metal/intermetallic  system  composed  of  titanium,  silicon,  and  tita¬ 
nium,  boron.  A  typical  shock  compaction  experiment  is  carried  out 
for  80%  porous,  stoichiometric  mixture  of  components,  with  the  ini¬ 
tial  mass  fractions  based  on  the  overall  equilibrium  products.  For 
the  Al,  CuO  system  the  stoichiometric  reaction  is  2  Al  +  3  CuO 
AI2O3  +  3  Cu.  The  reactive  material  sample  is  placed  in  a  striker  as¬ 
sembly,  and  compacted  by  the  action  of  two  metal  bars  that  are  shock 
loaded  by  the  firing  of  detonators  on  each  end.  A  sustained  heteroge¬ 
neous  front  was  found  to  propagate  at  an  average  speed  of  approxi¬ 
mately  6-20  cm/sec.  High  speed  microscopic  photography  was  used 
to  record  the  emitted  light  seen  through  a  small  observation  win¬ 
dow.  On  a  length  scale  of  10-20  p,m  one  observes  the  sudden  for¬ 
mation  and  disappearance  of  intense  spots  of  light,  corresponding  to 
intense  and  weak  chemical  reactions  recurring  within  a  time  interval 
of  approximately  100  |i  s.  Similar  results  were  observed  for  the  tita¬ 
nium  -  silicon  system  and  more  extensive  experiments  are  planned 
for  the  titanium  -  boron  system.  These  experiments  clearly  demon¬ 
strate  that  the  overall  combustion  process  is  highly  unsteady.  While 
the  lead  reactive  front  after  shock  compaction  is  observed  to  propa¬ 
gate  at  a  well-defined  average  velocity,  measurable  and  robust,  time- 
dependent  heterogeneous  reaction-diffusion  processes  occur  on  the 
micro-scale,  corresponding  to  the  initial  size  of  the  reactive  compo¬ 
nent  particles,  before  and  after  the  passage  of  the  lead  shock.  Since 
the  component  materials  and  their  reactants  are  very  hot  and  experi¬ 
ence  significant  thermal  expansion,  the  chemically  reacting  material 
experiences  a  distribution  of  local  flow  velocities  and  strain  rates,  pri¬ 
marily  at  the  material  interface  of  the  reacting  components. 


Fig.  1.  Schematic  of  the  condensed  phase  counterflow  diffusion  flame,  with  titanium 
entering  from  the  left  and  boron  from  the  right.  The  color  shades  clearly  delineate  the 
planar  mixing  region  and  reaction  zone. 

The  theoretical  modeling  approach  for  the  finely-space  laminate, 
or  regular  structured  materials,  can  be  summarily  described  as  an  ap¬ 
proach  that  lumps,  or  relies  on  cross-sectional  averages  for  all  trans¬ 
port  phenomena,  such  as  bulk  heat  transfer  and  diffusion,  and  for  all 
material  properties  of  the  laminate/arrays  [4].  These  reduced  mod¬ 
els  seem  to  appropriately  describe  observed  phenomena  that  de¬ 
pend  on  average  properties  of  the  system,  such  as  bulk  tempera¬ 
ture  and  reaction  extent,  or  the  self-propagation  speed  of  a  reactive 
front  propagating  along  the  axis  of  the  plane  of  the  plies.  They  do 
not  explicitly  describe  the  reaction-diffusion  phenomena  at  material 
interfaces,  and  cannot  delineate  separate  molecularly  distinct  reac¬ 
tants.  The  focus  in  this  work  is  exactly  on  the  processes  taking  place 
on  a  small-scale  of  the  initially-separated  component  materials  that 
comprise  the  mixtures,  that  are  not  necessarily  layered  or  structured. 
Our  approach  delineates  separate  molecularly  distinct  reactants  and 
employs  a  multicomponent,  thermodynamic  formulation  with  sepa¬ 
rated  reactants  and  products,  each  with  their  own  properties.  While 
we  employ  some  simplifications,  we  do  not  use  a  lumped,  averaged 
formulation  in  the  same  sense  of  the  reduced  models  found  in  the 
analysis  of  finely-spaced  laminates  or  arrays.  Fundamental  under¬ 
standing  of  these  local  events  will  serve  as  a  basis  for  future  model¬ 
ing  of  time-dependent  reaction  processes  in  both  classes  of  energetic 
materials,  agglomerated  composites  and  finely-spaced,  structured  or 
layered  composites. 

In  this  paper  we  examine  the  reaction-diffusion  processes  in  the 
mixing  region  of  two  opposed  streams  of  distinct  reactants,  that  stag¬ 
nate  on  the  stream  axis  as  shown  in  Fig.  1.  The  flow  due  to  mate¬ 
rial  deformation,  local  heating,  and  thermal  expansion  processes,  has 
been  presumably  established  at  the  micro-scale  by  prior  shock  com¬ 
paction  events  that  we  do  not  model;  such  considerations  were  given, 
for  example,  in  the  studies  reported  in  [7,9].  The  flow  of  opposed 
streams  is  characterized  by  a  single  strain  rate  -  the  ratio  of  the  av¬ 
erage  displacement  speed  to  the  length  scale  that  characterizes  the 
material  nonuniformity  -  or,  equivalently,  the  inverse  of  a  residence 
time  that  we  presume  given.  We  show  that  in  this  counterflow  con¬ 
figuration,  a  balance  between  advection,  diffusion  and  chemical  reac¬ 
tion  is  possible  under  steady  conditions  when  the  strain  rate  is  not  too 
large.  Otherwise,  the  only  possible  balance  is  that  of  a  nearly-frozen 
or  weak  burning  state  (not  discussed  in  this  paper).  In  the  larger  pic¬ 
ture,  we  envision  a  cooperative  time-dependent  mechanism  between 
local  reaction  sites  that  experience  intense  reaction  and  slow  reaction 
at  local  interfaces,  that  depends  on  the  local  strain  rate  distribution 
in  the  composite  material. 

Although  our  analysis  holds  for  any  three-components  sys¬ 
tem,  we  have  chosen  the  titanium  (Ti)  and  boron  (B)  system  as 
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a  representative  example  of  metal/intermetallic  reactive  material 
system.  Their  overall  reaction,  to  create  titanium  diboride  (TiB2) 
products, 

Ti  +  2B^TiB2,  (1) 

is  considered  under  nearly-isobaric  conditions,  such  that  the  pres¬ 
sure  is  represented  by  the  hydrostatic  pressure.  The  equation  of  state 
of  the  three  components  and  the  formulation  that  defines  the  equi¬ 
librium  equation  of  state  of  the  mixture  is  based  on  multicomponent 
thermodynamics  formulations  that  are  similar  to  those  used  in  the 
study  of  metallurgy  and  materials  [10],  As  such,  the  formulation  is 
based  on  Gibbs  thermodynamic  potentials  where  one  assumes  that 
at  each  point  of  the  condensed  phase  mixture,  there  is  a  single  stress 
state  and  temperature.  Each  isolated  component  is  assumed  to  have 
its  own  distinct  reference  density,  and  we  neglect  thermal  expansion 
in  the  components.  This  is  consistent  with  the  notion  that  the  change 
in  composition  due  to  reaction  is  much  larger  than  changes  due  to 
thermal  expansion.  As  a  result  the  mechanical  equation  of  state  for 
the  mixture  takes  a  simple  form  whereby  the  specific  volume  of  the 
mixture  is  simply  a  sum  of  the  intrinsic  densities  weighted  with  the 
mass  fraction  of  each  component.  This  form  of  the  mechanical  equa¬ 
tion  of  state  stands  in  contrast  with  that  for  a  mixture  of  reacting 
gases  that  is  a  relation  between  the  specific  volume,  pressure,  tem¬ 
perature  and  mass  fractions.  The  diffusion  model  for  the  components 
is  derived  from  an  effective  Fick  diffusion  formulation  as  described 
by  Curtiss  [11]  and  Curtiss  and  Bird  [12],  whereby  a  Maxwell-Stefan 
law  formulated  in  terms  of  binary  diffusivities  is  expressed  as  a  gen¬ 
eralized  Fick  diffusion  law  with  symmetric  diffusion  coefficients.  The 
resulting  diffusion  coefficients  used  in  our  model  are  then  chosen  to 
lead  consistent  results  with  experimental  data  reported  by  Trunov 
et  al.  [13], 

The  governing  equations  describing  the  steady  burning  of  initially 
separated  titanium  and  boron  in  a  counterflow  geometry  are  pre¬ 
sented  in  the  next  section.  These  equations  are  addressed  first  un¬ 
der  the  assumption  of  constant  mixture  density,  which  enables  the 
construction  of  an  analytical  solution  in  the  fast  chemistry  limit.  The 
expressions  for  flame  position  and  temperature  are  used  when  com¬ 
pared  to  experimental  data  to  estimate  binary  diffusion  coefficients 
that  are  not  well-known,  and  or  measured  under  restricted  condi¬ 
tions.  We  then  address  numerically  the  general  case,  of  finite  rate 
chemistry  and  when  the  density  of  the  mixture  varies  reflecting  the 
local  composition  of  the  mixtures.  In  particular,  we  derive  extinction 
conditions  that  differentiate  between  vigorous  and  weak  burning  be¬ 
tween  titanium  and  boron.  We  note  that  the  analytical  solution,  valid 
for  large  Damkohler  numbers,  is  used  to  validate  the  numerical  solu¬ 
tion  and  calibrate  the  corresponding  strain  rates  in  deriving  the  ex¬ 
tinction  conditions. 

2.  Formulation 

The  counterflow  geometry  under  consideration  is  shown  in  Fig.  1, 
where  far  to  the  left  (state  1)  there  is  only  titanium  and  far  to  the 
right  (state  2)  there  is  only  boron;  the  intrinsic  densities  are  de¬ 
noted  by  plo  and  p2o ,  respectively.  Under  steady  conditions,  the  ma¬ 
terial  deformation  may  be  described  by  a  velocity  field  v  of  the  form 
v  =  (u(x),  y  i>(x)).  This  “similarity  solution"  implies  that  the  pressure 
gradient  in  the  transverse  direction  y  is  necessarily  linear  and  admits 
planar  combustion  fronts  such  that  all  state  variables,  the  mass  frac¬ 
tions  Yj,  the  density  p,  and  temperature  T,  are  functions  of  x  alone. 
The  constituents  in  the  combustion  zone  include  titanium  of  mass 
fraction  Yj,  boron  of  mass  fraction  Y2,  and  titanium  diboride  prod¬ 
ucts  of  mass  fraction  Y3.  Overall  mass  conservation  implies  that 

Yi+Y2  +  Y3  =  1.  (2) 

The  conductivity  k  and  specific  heat  (at  constant  pressure)  cp  of 
the  mixture  (defined  as  mass-weighted  averages)  are,  in  general, 


functions  of  temperature  but  for  simplicity  they  will  be  taken  here 
as  constants.  Finally,  the  chemical  reaction  (1)  between  Ti  and  B  is 
assumed  to  proceed  at  a  rate 

a>  =  BYiY22  e~E/KT  (3) 

where  E  is  the  activation  energy,  Tl  is  the  universal  gas  constant  and  B 
is  an  appropriately  defined  pre-exponential  factor.  Different  reaction 
orders  could  be  considered  without  difficulty;  the  present  form  that 
simplifies  the  reaction  to  one-step  process  was  made  for  simplicity. 

2.3.  Conservation  equations 


The  governing  equations,  describing  conservation  of  mass,  mo¬ 
mentum,  and  energy  under  steady  conditions  simplify  to 


d  .  _ 

^(pu)  +  pv  =  0 

(4) 

dv  -2  _ 

puTx  +  pv  =C 

(5) 

dT  ,  d2j  n 

PUC^~kd^=(1C0 

(6) 

pi,<3x  +  5x(/0YlV,l)  =  -WlW 

(7) 

pLlW  +  ^pY2V2)  =  ~2W2CO  (») 

where  V,,  W;  stand  for  the  diffusion  velocity  and  molecular  weight  of 
species  i,  and  Q_  is  the  overall  heat  release.  As  noted  earlier,  the  pres¬ 
sure  gradient  in  the  transverse  y-direction  is  linear,  given  by  dp/dy  = 
-Cy,  where  C  is  a  constant  determined  by  the  far  field  boundary  con¬ 
ditions.  Specifically,  Eq.  (5)  implies  that  there  is  a  relation  between 
the  densities  and  strain  rates  at  the  far  ends,  such  that 

P\A  =  P2A=C-  (9) 

Hence  the  motion  of  the  reactants  impinging  on  each  other  is  char¬ 
acterized  by  a  single  strain  rate  e  (in  units  of  1/s)  which  we  choose 
as  e  =  €\/2\  the  factor  of  2  is  introduced  solely  to  facilitate  the  form 
of  the  analytical  asymptotic  solution  described  below.  The  axial  de¬ 
pendence  of  the  pressure  can  be  obtained  a-posteriori  by  solving  the 
x-component  of  the  momentum  equation  (not  written  above). 

Eqs.  (4)-(8)  must  be  supplemented  with  constitutive  relations  for 
the  diffusion  velocities  and  an  equation  of  state  for  the  mixture.  Start¬ 
ing  with  the  assumption  that  the  Gibbs  free  energy  of  each  com¬ 
ponent  of  the  mixture  can  be  summed,  the  Gibbs  potential  for  the 
mixture  which  is  a  function  of  pressure,  temperature,  and  mixture 
composition  is 

3 

g=X>(p.nVi.  (io) 

;=i 

where  the  energies  related  to  mixing  have  been  neglected.  The  spe¬ 
cific  volume  v  =  1/p  of  the  mixture  and  individual  components  are 
given  by  the  thermodynamic  relation: 


v  = 


dg 

dp 


and  Vj  = 


T.Y, 


d& 

dp 


(11) 


where  are  the  partial  volumes,  or  the  volume  of  the  components. 
This  leads  to 


,  =  J2V^P-  T)Yi 


(12) 


which  is  the  mechanical  equation  of  state  for  the  mixture.  (Note:  The 
symbol  v  represents  the  specific  volume  and  should  not  be  confused 
with  v  which  represents  the  y-component  of  velocity.)  If  we  further 
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assume  that  the  volume  change  due  to  pressure  variations  is  small, 
and  neglect  the  effect  of  temperature,  the  partial  volumes  u,  can  be 
approximated  by  their  reference  values  0,Q.  The  mechanical  equation 
of  state  (12),  when  expressed  in  terms  of  the  densities,  then  becomes 


are  concentration  dependent  and  may  not  necessarily  be  all  positive. 
They  must,  however,  satisfy  the  constraints  [12] 

=  for  all  i,  j,  ^BjjY;  =  0  for  all  j. 


p~'  =  EYi  A;1 


(13) 


where  pio  is  the  intrinsic  density,  and  p^1  the  intrinsic  specific  vol¬ 
ume,  of  species  i.  Using  (2)  the  equation  of  state  simplifies  to 


With  the  flow  specified  by  the  strain  rate  e  and  the  material  prop¬ 
erties  far  to  the  left  (denoted  by  subscript  1)  and  far  to  the  right 
(denoted  by  subscript  2)  assumed  uniform,  the  boundary  conditions 
are: 

du/dx - 2e  as  |x|  -»  oo,  (15) 

P  =  Pi0,  Yj  =  1,  Y2  =  0.  T  =  Too  as  -oo,  (16) 

P  =  P20.  Y\  =  0,  Y2  =  1,  T  =  Too  as  x->  +oo.  (17) 

Note  that  there  is  no  need  to  specify  a  condition  for  v,  because  it  is 

obtained  from 

~v=~-pTx{pu)  (18) 

by  differentiation. 

When  density  variations  are  small,  the  equation  of  state  (14)  can 
be  effectively  replaced  by  p  =  constant,  and  the  velocity  field  every¬ 
where  is  given  by 

u  =  -2ex.  v=2e.  (19) 

The  problem  reduces  to  the  reaction-diffusion  system  (6)-(8).  The 
constant-density  approximation  will  be  used  for  simplicity  in  the 
asymptotic  description  described  below.  In  general,  variations  in 
the  density  modify  the  overall  velocity  field  which,  in  turn,  affects 
the  combustion  field.  In  Section  4,  numerical  computations  are  car¬ 
ried  out  in  order  to  assess  the  importance  of  density  variations  in 
condensed-phase  combustion. 


P  P30  \Pl0  P3„ 


2.2.  Diffusion 

The  most  common  expressions  used  for  multi-component  diffu¬ 
sion  are  the  Maxwell-Stefan  (MS)  relations 

Y.Y. 

(20) 

i  Ai 

with  the  summation  taken  over  all  species  present  [14];  here  X,  is 
the  molar  fraction,  Vj  is  the  diffusion  velocity  vector  of  species  i,  and 
@y  =  3>ji  is  the  binary  diffusivity  of  a  pair  of  species  (i,  j).  Although 
this  relation  was  derived  for  a  dilute  ideal  gas  mixture,  it  has  been 
often  applied  to  condensed  phase  media  [15], 

The  use  of  the  Maxwell-Stefan  relations  is  quite  complicated  be¬ 
cause  the  diffusion  velocities  Vj  are  not  expressed  explicitly  in  terms 
of  the  concentration  gradients.  A  common  practice  is  to  use  the  gen¬ 
eralized  Fick  equations 

Vi^ByVX,-  (21) 

j 

with  coefficients  By,  referred  to  as  Fick  diffusivities ,  that  are  related 
to  the  binary  diffusivities  @y,  but  unlike  the  binary  diffusivities  they 


The  expressions  relating  Fick  and  binary  diffusivities  for  a  ternary 
mixture  [  16]  are  listed  in  the  Appendix.  Converting  Eq.  (21 )  from  mole 
to  mass  fraction  yields: 

YiV i  =  a.-VY,  +  5iVY2,  i=  1,2  (22) 

where 


way2  +  y3)(~y2w3  -  w2+Y2W2)@ng>u  +  Y,Y2w2(Wi  -w3)@23®n 

YlW2W3@23  +  Y>W,W3®13  +YsW2W,®12 

W3 Y,  (W2  -Y2W2  +  Y>  W,  )@23  @, 3 
Y,  W2  W3  @23  +  Y2  W,  W3  @13  +  y3  W2W,  3>n 

YiW]  (y2  +  Y3)(W2  -  W3)@13@12  +  YiWif-W]  +  YiWi  -WjY,)^®^ 
Y]  W2W3@23  +  y2W,  W3@,  3  +  YsM^W,  @,2 
Y,  W3  (— Y,  VV2  —  W,  +  Y,  W, )  @23  @13 
+  Yi  W2W3  @23  +  YAV,  W3  @13  +  Y3W2Wi  @12 
Y2W,  (-VV2  +  Y2W2  -  Y2W3)@13@12  +  Y2W2(Yi  +  Y3)(W,  -  W3)@23@i2 

y3  W2W3  @23  +  y2w,  w3  @13  +  y3  vv2  w,  @,  2 

Y2W2(-Y2W2 +Y2W^  +  W2)@23@13 
“  Y, W2W3@ 23  +  YiW] W3@13  +  Y3W2W1  @,2 

Y2Y1W3  (— W3  +  W2)@i3@i2  +  Vy2(Yi  +  Y3)(— Wi  +Y,W,  -  Y1W3)@23@12 
Y,  W2 W3  @23  +  Y>  W,  VV3  @33  +  Y3 W2W\  @, 2 
Y2W3(-Y1W2  -  W,  +Y,W1)@23@13 
+  Y,  W2W3  @23  +  Y2Wi  W3@13  +  YjWjW,  @32  ' 


The  species  Eqs.  (7)-(8)  can  then  be  written  as 
dY,  d 


1  S 


K 


dY1  ,  dY2\ 
ai^TT  +  bi^r) 


dx 


dx 


=  -  W1  co, 


(23) 


pu 


dY2 

dx 


d 

dx 


p[a2 


dY1 

dx 


dY2 
'  dx 


=  -2W2  co. 


(24) 


A  simplification  that  can  be  used  for  analytical  convenience  results 
from  assuming  equal  molecular  weights  =W2  =  W3,  then 


@12  +  (@23  -  @12)^1 

d  1  —  "4  - 

(@23  -  @12) Yi  +  (@13  -  @i2)Y2  +  @i2 

^  _  _ @23  (@12  ~  @13)^1 _ 

1  “  (@23  -  @12)^1  +  (@13  -  @12)^2  +  @12 

a  =  _ @13  (@12  ~  @23)^2 _ 

(@23  -  @12)^1  +  (@13  -  @12)^2  +  @12 
,  @12  +  (@13  -  @12)^2 

m  =  — - 

(@23  -  @12) Yi  +  (@13  -  @12)Y2  +  @12 


3.  Asymptotic  solution  -  the  Burke-Schumann  limit 

We  first  present  analytical  results  in  the  limit  of  infinitely  fast 
chemical  reaction,  known  in  the  literature  as  the  Burke-Schumann 
limit.  The  solutions  obtained  in  this  limit  provide  a  simple  illustra¬ 
tion  of  the  flame  structure.  An  appropriate  time  scale  may  be  defined 
as  the  ratio  of  a  representative  length  associated  with  the  distance 
between  particles  and  a  characteristic  microscale  velocity  associated 
with  material  deformation.  The  inverse  of  this  time  is  the  charac¬ 
teristic  strain  rate,  and  the  fast  chemistry  limit  corresponds  to  weak 
strain  rates.  The  corresponding  Damkohler  number,  or  the  ratio  of  the 
flow-to-chemistry  time  scales,  is  therefore  large.  Although  asymp¬ 
totic  methods  that  accounts  for  finite  rate  chemistry  and  thus  span 
a  wider  range  of  strain  rates  (up  to  flame  extinction)  are  available  for 
the  related  gaseous  problem  [17,18],  their  extension  to  energetic  ma¬ 
terials  is  nontrivial  and  will  be  discussed  in  a  future  publication.  Here 
we  rely  on  numerical  methods  to  examine  the  dependence  of  the  so¬ 
lution  on  the  strain  rate  for  steady  combustion. 
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In  the  fast  chemistry  limit  the  chemical  reaction  occurs  along  a 
sheet,  at  x  =  Xf  say,  where  the  two  reactants  are  in  contact.  Else¬ 
where,  the  chemical  reaction  is  negligibly  small  and  we  are  left  solv¬ 
ing  the  energy  and  species  equations  on  either  side  of  the  sheet,  with 
co  =  0.  The  flame  sheet  separates  a  region  where  there  is  only  tita¬ 
nium  ( x  <  Xf),  from  a  region  where  there  is  only  boron  (x  >  Xf).  Since 
Y2  =  0  for  x  <  Xf,  we  find  that 


ai  =  -®ri, 

a2  =  0, 


^23(^12  -  @u)Yi 

5*23^1  +  ^12(1  —  Yi) 

^23^12 

^23^1  +  @120  —  ¥1) 


and,  since  Y,  =  0  for  x  >  Xf,  we  find  that 

^13^12 


a  1  =  - 

a2 


^13^2  +  ®12(1  -  V2) 
^13(^12  -  ^23)^2 
^13^2  +  ^12(1  “  Y2) 


hi  =0. 

b2  =  —®22. 


We  note  parenthetically  that  the  simplification  of  the  coefficients  a,, 
fa,  applies  even  for  unequal  molecular  weights.  All  variables  must  be 
continuous  at  the  flame  sheet,  but  their  derivatives  are  not;  the  mass 
and  energy  fluxes  must  satisfy  the  relations 


a 

dT 

1 

dY,  , 

dY2 

1 

Q./  cp 

dx 

chTC  +  b' 

dx 

~  2SN2 

(25) 


where  the  square  brackets  [■]  denote  the  jump,  namely  the  difference 
between  the  values  atxj:  and  xj,  and  a  =  k/pcp  is  the  thermal  diffu- 
sivity  of  the  mixture.  The  conditions  (25)  imply  that  the  fluxes  of  ti¬ 
tanium  and  boron  towards  the  flame  sheet  are  in  stoichiometric  pro¬ 
portions,  and  they  specify  the  proportion  of  heat  from  the  total  heat 
released  conducted  to  one  or  the  other  side  of  the  sheet. 

For  simplicity,  we  have  also  adopted  in  this  section  the  constant 
density  approximation.  The  mathematical  problem  on  either  side  of 
the  flame  sheet  then  consists  of 


2ex-j — ha-r^-=0  forx^gXf  (26) 

ax  ax2  1 

dY  d2Y 

2ex-j —  +  ^13  —7-^-  =  0  for  x  <  xf  (27) 

ax  ax 2  J 

Yi  =  0  for  x  >  Xf  (28) 

Y2  =  0  for  x  <  Xf  (29) 

dY  d2Y 

2ex-j-^  +  ^23-j-^  =  0  forx>x^  (30) 


together  with  the  boundary  conditions  (16)  and  (17),  with  p  assumed 
constant,  and  the  jump  relations 


[T'l  =  [Vi]  =  [y2]  =  0 


(31) 


a 

dT 

®13 

dY, 

®23 

Q/Cp 

dx 

W, 

dx 

2  W2 

dx 

across  x  =  Xf.  The  solution  is  readily  obtained  as 


T  = 


Too  +  (J/  —  Too) 


Too  +  ( Tf  —  Too) 


1  +  erf(Ve/Q!  x) 
1  +erf (ye /a  xf ) 

1  -  erf [sI^Iol  x) 

1  -  er f^e/axy) 


X  <  Xf 


X  >  Xf 


(32) 


Fig.  2.  Variation  of  the  flame  sheet  position  with  ^23/^13,  based  on  the  asymptotic 
solution. 


Y, 


Y2 


1  +  erf(Ve/®i3  x) 
1  +erf(Je/@uxf) 


1  -  erf(yg/g>23x) 
1  -  erf (ofcWnXf) 


X  <  Xf 
X  >  Xf 
X  <  Xf 
X  >  Xf 


The  position  Xf  of  the  flame  sheet  and  the  adiabatic  flame  tempera¬ 
ture  Tf,  defined  as  the  value  of  T  at  the  flame  sheet,  satisfy 

1  +erf(Ve/^13x/)  =  e£*r/a)23 

1  -  erf {ofcJW^yXf)  V Y  ®23  gCxj/^n 


1  Q/Cp  fa  1-erf 2(V^xf)  e^'a 
f  00  2  W]  V  “  1  +  erf (s/e/9ri  X/)  e«J/®,3 


(34) 


where  v  =  2\A/2/W-[  is  a  mass-weighted  stoichiometric  ratio.  The  po¬ 
sition  Xf  is  determined  from  the  transcendental  Eq.  (33)  using  an  it¬ 
erative  process.  We  note  that  for  a  given  strain  rate  e  the  position  Xf 
depends  only  on  the  binary  diffusivities  Ti-TiB2  and  B-TiB2,  and  is  in¬ 
dependent  of  the  diffusivity  of  Ti-B  because  there  is  no  boron  in  the 
titanium  region  and  vice-versa.  Once  Xf  is  determined,  the  flame  tem¬ 
perature  can  be  calculated  from  (34)  by  direct  evaluation.  Evidently, 
the  latter  depends  on  the  heat  released  Qand  thermal  diffusivity  a. 

For  equal  diffusivities  ®\3  =  ®23  =  Eq.  (33)  reduces  to 

erf (y/e/9Xf)  = 

For  the  titanium-boron  reaction  the  mass-weighted  stoichiometric 
ratio  v  ~  0.45,  implying  that  Xf  -0.41  y^/e  and  the  flame  sheet 
lies  on  the  titanium  side  of  the  stagnation  plane.  If  the  Lewis  number 
is  assumed  equal  to  one,  i.e.,  ®  =  a,  the  flame  temperature  is  given 
by 


Tf  —  Too  + 


Q/CpVVj 
1  +  v  ' 


In  the  absence  of  differential  (unity  Lewis  number)  and  preferential 
(unequal  mass  diffusivities)  diffusion,  the  flame  temperature  results 
from  a  simple  energy  balance. 

Figure  2  displays  the  dependence  of  the  scaled  flame  sheet  loca¬ 
tion  on  the  ratio  ®i23l®i\3,  for  a  fixed  v.  This  ratio  is  typically  larger 
than  one,  since  boron  atoms  have  an  effectively  smaller  atomic  ra¬ 
dius  and  hence  diffuse  more  readily  through  titanium-diboride  than 
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does  the  titanium  atom  [19,20].  We  see  that  the  flame  sheet  gener¬ 
ally  resides  on  the  titanium  side  of  the  stagnation  plane  and  moves 
further  to  the  left  when  the  diffusivity  of  boron  into  the  titanium- 
diboride  products  increases  relative  to  the  diffusivity  of  titanium  into 
the  titanium-diboride  products. 

The  only  explicit  results  presented  here  from  the  asymptotic  solu¬ 
tion  are  the  expressions  for  the  flame  position  and  temperature.  Pre¬ 
dictions  based  on  these  expressions  are  used  in  conjunction  with  ex¬ 
perimental  data  to  extract  kinetic  parameters  that  are  not  well  known 
or  were  measured  under  limited  conditions.  Moreover,  the  asymp¬ 
totic  solution  is  used  to  validate  the  numerical  solution  of  the  gov¬ 
erning  equations  for  large  values  of  the  Damkohler  number  (or  small 
strain  rates)  and  to  properly  calibrate  the  strain  rate,  in  order  to  derive 
extinction  conditions. 

4.  Numerical  solution  -  finite-rate  chemistry 


Table  1 

Property  values  of  representative  physical  parameters  used  in  the  computa¬ 
tions. 


Property 

Symbol 

Value 

Reference 

Heat  release 

a 

—323.8  kj/mol 

[23] 

Molar  mass  ofTi 

Wl 

47.87  g/mol 

Molar  mass  of  B 

w2 

10.81  g/mol 

Molar  mass  ofTiB2 

w3 

69.85  g/mol 

Averaged  molar  mass 

wc 

50  g/mol 

Intrinsic  density  of  Ti 

Ph 

4.5  g/cc 

Intrinsic  density  of  B 

P2o 

2.34  g/cc 

Intrinsic  density  of  TiB2 

P3o 

4.52  g/cc 

Averaged  density 

Pco 

3.8  g/cc 

Heat  capacity 

CP 

900  J/(kg  I<) 

[23] 

Thermal  conductivity 

k 

36  W/(m  K) 

[24-26] 

Pre  exponential  factor 

B 

7.6el6  mol/(m3  s) 

[27] 

Activation  energy 

Ea 

318  kj/mol 

[27] 

Binary  diffusivity  of  Ti-B 

®12 

0.2  m2/s 

[28] 

To  examine  the  effects  of  finite-rate  chemistry,  the  boundary 
value  problem  consisting  of  Eqs.  (4)-(6)  and  (23)-(24)  and  bound¬ 
ary  conditions  (15)-(17)mustbe  addressed.  The  numerical  procedure 
is  described  next,  followed  by  results  pertaining  to  titanium-boron 
combustion. 


4.3.  Numerical  procedure 


Steady  solutions  of  the  aforementioned  boundary  value  problem 
are  obtained  numerically  as  the  long  time  behavior  of  the  time- 
dependent  equations 


dv 


dv 


pdt+puTx 


■  pi)2  =  C 


(35) 


9Y, 

9Y, 

9 

r  / 

9Y,  . 

3Y2\ 

=  -Wj« 

pHF 

+ 

pUJt  + 

dx 

P/Qi 

-rT~  +  b 

dx 

dx  ). 

9Y2 

9Y2 

9 

r  / 

9Y,  , 

9Y2\7 

=  —2 W2  co 

P^t 

+ 

pU~d^  + 

dx 

P[a  2 

-TC—  +  i>2 

dx 

dx  ). 

9T 

dT\ 

92T 

Qoo. 

Pcp[ 

~dt 

+  U  9x) 

-  k 

9X2  “ 

with  v  obtained  from 

pv=~±(pu) 


(36) 

(37) 

(38) 


(39) 


and  p  from  (14).  initially,  these  equations  were  solved  using  an  ex¬ 
plicit  time  marching  method  until  the  solution  converges  to  its  equi¬ 
librium  state.  The  integration  in  time  starts  with  an  initial  guess, 
taken  here  as  the  asymptotic  solution  discussed  in  the  previous  sec¬ 
tion.  A  fourth  order  approximation  was  used  to  compute  the  first  or¬ 
der  and  second  order  space  derivatives  and  a  fourth  order  Runge- 
Kutta  (RK4)  method  was  used  for  time  stepping.  The  extent  of  the 
numerical  domain  depends  on  the  strain  rate  value,  with  lower  strain 
rates  requiring  a  larger  domain.  This,  however,  can  be  overcome  by 
normalizing  x  with  the  thermal  diffusion  length  ld  =  y/a/e.  as  is  also 
evident  from  the  analytical  form  of  the  asymptotic  solution.  Due  to 
large  stiffness  arising  from  the  Arrhenius  exponential  in  the  reaction 
rate  term,  we  found  that  the  time  step,  in  general,  could  not  exceed 
10~9  and  to  properly  describe  the  solution  at  low  strain  rates  where 
the  reaction  zone  becomes  extremely  thin,  a  fine  grid  is  also  required. 
As  a  result,  convergence  to  the  equilibrium  state  was  very  slow  even 
after  parallelizing  the  numerical  code.  Determining  the  solution  over 
a  wide  range  of  strain  rate  conditions  requires  a  faster  converging 
algorithm. 

To  overcome  the  computational  stiffness,  an  implicit  time  re¬ 
laxation  method  was  implemented.  First  and  second  order  spa¬ 
tial  derivatives  are  approximated  by  second  order  finite  difference 


schemes  on  a  uniform  grid,  and  for  time  stepping  we  use  a  back¬ 
ward  Euler  method,  such  that  a  generic  equation  of  the  form  dcp/dt  = 
f(t,<p)  is  approximated  by  0n+1  =  0n  +  At/(tn+1,(/>"+1),  where  n 
denotes  the  time  step  and  At  the  time  increment.  A  damped  Newton- 
Raphson  solver  is  then  used  for  the  solution  of  the  resulting  nonlin¬ 
ear  system  at  each  time  step.  The  stiffness  in  the  governing  equa¬ 
tions  lead  to  strong  fluctuations,  and  consequently  sharp  gradients 
in  the  Jacobian  matrix  that  may  cause  the  numerical  solution  to  di¬ 
verge.  We  have  therefore  implemented  a  relaxation  parameter  dur¬ 
ing  each  iterative  update  to  limit  the  fluctuations;  the  value  0.1  was 
found  appropriate  for  our  calculations.  The  PetSC  sparse  solver  [21] 
was  used  for  the  system  that  arises  during  the  Newton-Raphson  iter¬ 
ation.  This  approach  has  significantly  reduced  the  time  step  relative 
to  the  explicit  scheme  from  10~9  to  10  2 ,  for  the  same  spatial  grid 
distribution. 

In  order  to  generate  solutions  for  increasing  values  of  the  strain 
rate  e  and  draw  response  curves  of  quantities  of  interest  (e.g.,  flame 
temperature,  mass  fraction  of  unconsumed  reactants,  etc.)  as  a  func¬ 
tion  of  6,  we  start  with  a  small  strain  rate  value  of  e  =  0.01  s"1  using 
the  asymptotic  solution  as  an  initial  guess,  and  advance  in  time  until 
the  incremental  changes  in  the  solution  at  all  points  in  the  domain 
of  integration  are  less  than  a  tolerance  error,  here  taken  as  10-3.  Con¬ 
vergence  is  guaranteed  only  when  the  long-time  equilibrium  solution 
is  stable.  This  approach,  therefore,  works  well  for  small-to-moderate 
values  of  e,  but  fails  at  larger  values  when  the  solution  of  the  bound¬ 
ary  value  problem  becomes  multi-valued,  and  the  response  curve  de¬ 
velops  a  turning  point  with  stable  and  unstable  branches.  To  con¬ 
struct  all  steady  solutions,  stable  and  unstable,  particularly  near  the 
turning  point,  we  use  the  continuation  approach  proposed  by  Kur- 
dyumov  and  Matalon  [22],  The  time-dependent  Eqs.  (35)-(39)  were 
solved  with  an  additional  constraint  that  the  temperature  remains 
constant  at  some  reference  point,  say  T(x*)  =  T*.  The  constraint  is 
used  to  iterate  on  the  value  of  e  until  both,  the  strain  rate  and  the 
space  distribution  of  solution,  do  not  vary  significantly  from  one  time 
step  to  the  next.  By  selecting  Vs  judiciously,  this  procedure  always 
converges  allowing  to  generate  the  entire  response  curve  through  a 
turning  point.  When  computing  steady  states  that  are  expected  to  be 
unstable,  the  time  steps  must  be  adaptively  changed  while  keeping 
the  damped  newton  iterations.  It  is  crucial  to  start  with  small  time 
steps,  10~6  say,  and  gradually  increase  At  to  10-2.  Finally,  we  note 
that  the  time  t  is  used  in  this  context  only  as  an  iterative  parameter. 

4.2.  Physical  parameters 

Table  1  lists  representative  values  of  physical  parameters  based  on 
a  literature  survey;  some  values  required  by  the  model  are  easier  to 
estimate  than  others  and  are  obtained  from  standard  thermal  prop¬ 
erty  measurements.  The  table  lists  values  for  the  constant  pressure 
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(a)  (b) 

Fig.  3.  Variation  of  the  adiabatic  flame  temperature  7)  with  (a)  variations  in  ®i3,  forgiven  ®2  3  =  5.0  ■  10  6  m2/s,  and  (b)  variations  in  ®23,  forgiven  ®]3  =  5.0  10  7  m2/s,  relative 
to  base-line  values  marked  by  a  dashed  line. 


(a)  Temperature  profile 


(b)  Mass  fractions  profiles 


Fig.  4.  Comparison  between  the  numerical  solution  and  asymptotic  solutions  for  a  small  strain  rate  value  e  =  0.01  s_1 . 


heat  capacity,  cp,  and  thermal  conductivity,  k,  that  represent  aver¬ 
aged  values  for  the  mixture.  However,  the  determination  of  the  bi¬ 
nary  mass  diffusivities  is  more  problematic.  Experimental  values  for 
3  and  ^23  have  been  reported  [19,20],  but  not  at  conditions  that 
are  present  in  the  reaction  zone.  Better  estimates  could  be  obtained 


from  our  derived  asymptotic  expressions,  when  used  in  conjunction 
with  experimentally  measured  flame  temperatures. 

Trunov  et  al.  [13]  measured  adiabatic  flame  temperatures  for 
the  titanium-boron  reaction  of  approximately  2400-3300  K.  One 
expects,  similar  to  the  values  measured  experimentally  at  low 
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(a)  Flame  position 


(b)  Flame  temperature 


Fig.  5.  Response  curves  of  flame  position  and  flame  temperature  versus  strain  rate. 
The  two  curves  (red/blue)  correspond  to  the  complete  and  simplified  diffusion  formu¬ 
lations.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is 
referred  to  the  web  version  of  this  article.) 

temperatures  [20],  that  the  binary  diffusivities  of  both,  boron  and 
titanium  into  titanium-diboride,  are  significantly  smaller  than  the 
mixture  thermal  diffusivity,  and  that  @23  »  £^13  •  Taking  ^23/^13  = 
10  as  a  base-line  model  value,  and  the  estimate  Tf  =  3000  K  for 
the  flame-sheet  temperature,  our  formulas  (33)  and  (34)  yield  ®13 
5  x  10-7  m2/s  and  S>22  ~  5  x  10~6  m2/s.  Figure  3  shows  changes  in 
flame  temperature  as  a  result  of  variations  in  the  binary  diffusivi¬ 
ties  ®13  and  £%3  relative  to  the  base-line  values,  for  the  fixed  ratio 
of  ^23/^13- 

4.3.  Low  strain  rates 

We  start  by  presenting  results  pertaining  to  low  strain  rates, 
where  the  solution  can  be  compared  directly  to  the  asymptotic  so¬ 
lution  discussed  above.  Since  the  latter  was  obtained  under  the  con¬ 
stant  density  assumption,  the  average  value  p  =  3.8  g/cc  was  se¬ 
lected,  abandoning  the  equation  of  state  (14).  The  expressions  for  the 
diffusion  coefficients  were  used  without  resorting  to  the  approxima¬ 
tion  of  equal  molecular  weights,  since  the  reduced  form  of  these  re¬ 
lations  led  to  nearly  identical  results. 

Figure  4  (a)  and  (b)  shows  a  comparison  of  the  temperature 
and  mass  fraction  profiles  between  the  computed  solution  for  e  = 
0.01  s_1  and  the  corresponding  asymptotic  expressions.  The  spa¬ 
tial  coordinate  is  normalized  with  the  thermal  diffusion  length  ld  = 
~Ja/e  =  3.2410~4  cm.  For  such  low  values  of  strain  rate  combustion 


(a)  Unconsumed  titanium  mass  fraction 


(b)  Unconsumed  boron  mass  fraction 


Fig.  6.  The  extent  of  unconsumed  reactants  leaking  through  the  reaction  zone.  The  two 
curves  (red/blue)  correspond  to  the  complete  and  simplified  diffusion  formulations. 
(For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred 
to  the  web  version  of  this  article.) 

is  nearly  complete,  with  both  reactants  consumed  in  a  very  thin  reac¬ 
tion  zone.  There  is  excellent  agreement  between  the  computed  and 
asymptotic  profiles  except  for  very  small  changes  near  the  reaction 
zone,  as  shown  in  the  figure  inserts.  Evidently,  for  finite  e  however 
small,  the  reaction  zone  has  a  finite  thickness,  which  is  adequately 
captured  by  the  numerical  solution  as  shown  in  Fig.  4(c).  The  excel¬ 
lent  agreement  between  the  numerical  and  asymptotic  solutions  also 
serves  as  a  validation  of  our  numerical  methodology  that  properly 
resolves  the  stiffness  in  the  governing  equations  arising  from  the  ex¬ 
ponential  Arrhenius  term. 

4.4.  Moderate  strain  rates  -  constant  density 

Next  we  consider  the  entire  response  of  the  flame  to  increasing 
strain  rates,  from  complete  combustion  corresponding  to  low  strain 
rates  to  flame  extinction  occurring  at  significantly  higher  values  of  e. 

We  first  examine  the  solution  using  two  diffusion  formulations: 
the  complete  expressions  for  the  diffusion  coefficients  a,-,  and  the 
simplified  form  resulting  from  equal  molecular  weights.  We  note 
that  the  general  diffusion  expressions  add  nonlinearity  to  an  already 
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Fig.  7.  Profiles  of  the  transverse  velocity  v  across  the  combustion  zone  for  two  strain 
rate  values. 
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Fig.  8.  Density  profiles  across  the  combustion  zone  for  two  strain  rate  values. 


stiff  problem.  Therefore,  to  facilitate  the  computations  we  have  re¬ 
tained  the  constant-density  approximation  in  the  results  presented 
in  this  subsection,  which  effectively  decouples  the  flow  and  combus¬ 
tion  fields;  variable  density  solutions  will  be  presented  in  the  next 
subsection.  Figures  5  and  6  show  response  curves  of  various  flame 
characteristics  to  the  imposed  strain  rate  e,  based  on  the  two  dif¬ 
fusion  formulations.  These  include  the  flame  position  Xf  (defined  as 
the  location  where  the  temperature  reaches  its  maximum  value),  the 
flame  temperature  Tf  (the  value  of  T  at  Xf),  and  the  mass  fractions 
of  unconsumed  titanium  and  boron  (the  values  of  Y ^  and  Y2  evalu¬ 
ated  at  xf).  We  note  that  the  precise  evaluation  of  the  flame  position, 
and  the  temperature  and  mass  fractions  evaluated  at  this  location, 
depend  on  the  step  size  used  in  the  computations  and  on  the  value  of 
the  thermal  diffusion  length  ld,  which  decreases  from  3.24  x  10~5  cm 
to  1.3  x  10~6  cm  as  the  strain  rate  increases  from  e  =  0.01  s_1  to 
e  =  6.5  s_1  (which  is  very  near  extinction).  A  Matlab  Curve  Fit  [29] 
tool  was  used  to  fit  the  numerical  data  to  a  smoothing  spline,  which 
was  then  used  for  the  determination  of  Xf  and  for  evaluating  all  prop¬ 
erties  at  this  location.  The  same  tool  was  used  to  fit  smoothing  spline 
curves  presented  in  all  the  figures,  to  the  discrete  data  obtained 
numerically. 

From  the  response  curves  of  Figs.  5  and  6  the  following  physical 
picture  emerges.  Since  at  low  strain  rates  the  chemical  reaction  time 
is  much  shorter  than  the  flow  time,  the  reaction  proceeds  immedi¬ 
ately  as  titanium  and  boron  get  in  contact.  The  reaction  occurs  in  a 
very  thin  zone  (or  a  sheet)  where  both  reactants  are  completely  con¬ 
sumed.  The  flame  lies  on  the  titanium  side  (x  <  0)  and  the  flame  tem¬ 
perature  reaches  its  maximum  value.  Upon  increasing  the  strain  rate, 
the  flow  time  relative  to  the  chemical  reaction  time  is  shortened  and 
a  fraction  of  titanium  and  boron  escape  mixing  and  leaks  through 
the  reaction  zone.  As  a  result  of  incomplete  combustion,  the  flame 
moves  towards  the  boron  side  and  the  flame  temperature  drops.  The 
relatively  larger  leakage  of  titanium  as  opposed  to  boron  stems  from 
the  fact  that  gq3  <g)  S>22,  which  implies  much  larger  fluxes  of  boron 
towards  the  reaction  zone  and  consequently  a  more  complete  con¬ 
sumption  of  boron.  This  trend  continues  until  the  fraction  of  uncon¬ 
sumed  reactants  exceeds  a  critical  threshold  and  the  flame  temper¬ 
ature  drops  significantly  from  its  adiabatic  value.  For  larger  values 
of  the  strain  rate  steady  burning  is  no  longer  possible.  The  critical 
state,  represented  by  the  turning  point  on  the  response  curves  and 
corresponding  to  e  ^  6.586  sec^1,  identifies  conditions  associated 
with  flame  extinction.  The  lower  branch  on  the  temperature  response 
curve  of  Fig.  5  and  the  upper  branches  in  Fig.  6  are  typically  unstable 
and  therefore  physically  inaccessible. 


(a)  flame  psotition 


(b)  flame  temperature 


Fig.  9.  Comparison  of  flame  position  and  temperature  between  constant  and  variable 
density  conditions. 

The  two  diffusion  formulations  lead  to  identical  results  at  low 
strain  rates  and  predict  the  exact  same  extinction  strain  rates.  There 
are  small,  insignificant  differences  in  flame  temperature  at  and  near 
extinction,  which  can  be  traced  to  the  slight  difference  in  flame  loca¬ 
tion.  Being  influenced  primarily  by  the  binary  diffusivity  5?13  <sc  ^23. 
the  simplified  diffusion  formulation  predicts  a  flame  position  that  is 
slightly  tilted  towards  the  titanium  (x  <  0)  side.  Due  to  the  negligible 


S.P.  Koundinyan  et  al.  /  Combustion  and  Flame  162  (201 5)  4486-449 6 


4495 


(a)  Unconsumed  titanium  mass  fraction 


(b)  Unconsumed  boron  mass  fraction 

Fig.  10.  Yi  response  curve  comparison  between  constant  and  variable  density 
approximations. 

difference  between  the  two  formulations,  the  simplified  diffusion  for¬ 
mulation  will  be  used  in  the  following  section  in  order  to  save  com¬ 
putational  time. 

4.5.  Moderate  strain  rates  -  variable  density 

When  variations  in  density  are  accounted  for,  the  flow  field  no 
longer  satisfies  (19)  and  must  be  obtained  by  solving  the  momentum 
Eq.  (5)  with  p  given  by  (14).  The  boundary  conditions  as  x  ->  ±oo 
imply  that  v  tends  to  a  constant  on  either  end,  with  the  asymptotes 
subject  to  the  constraint  (9).  Figure  7  shows  profiles  of  v  computed  for 
two  values  of  e;  the  low  strain  rate  value  corresponding  to  conditions 
close  to  the  Burke-Schumann  solution  and  the  larger  strain  rate  value 
corresponding  to  near-extinction  conditions.  Variations  in  v  across 
the  field  are  due  to  density  changes,  which  are  more  pronounced  in 
the  combustion  region  -0.5  <  x/ld  g  0.5. 

Density  profiles  are  shown  in  Fig.  8  for  the  same  two  values  of 
e.  Generally,  the  density  would  vary  from  plo  =  4.5  g/cc  on  the  tita¬ 
nium  side  to  p2o  =  2.34  g/cc  on  the  boron  side.  Since  to  the  density 
of  TiB2  is  slightly  higher  than  that  of  Ti,  there  is  a  noticeable  jump  in 
density  near  the  reaction  zone  at  low  strain  rates  (see  the  figure  in¬ 
set),  where  the  reaction  is  more  vigorous  and  the  production  of  TiB2 
higher.  This  increase  in  density  diminishes  at  higher  strain  rates  due 
to  lower  TiB2  production  and  increase  of  titanium  leakage  through 
the  reaction  zone. 

Response  curves  of  flame  position  and  temperature  vs  strain  rates, 
for  constant  and  variable  density  conditions,  are  shown  in  Fig.  9.  The 


flame  temperature  Tf  and  flame  position  Xf  are  consistently  higher  for 
the  constant  density  case  than  for  the  variable  density  case,  for  all  val¬ 
ues  of  e.  The  difference  is  very  small  and  is  of  no  physical  significance, 
being  primarily  due  to  the  selected  mean  density  value  adopted  in  the 
constant  density  formulation.  The  small  difference  between  the  two 
solutions  indicates  that  the  composition  distribution  has  little  effect 
on  the  overall  density,  for  practically  all  strain  rate  values.  The  flame 
temperature  at  extinction  is  approximately  2500  K,  corresponding  to 
a  drop  in  approximately  500  K  from  the  adiabatic  flame  temperature 
Tf  =  2950  K  (or  3000 1<  for  the  constant  density  case).  The  flame  posi¬ 
tion  in  the  Burke-Schumann  limit  is  determined  from  the  condition 
that  the  two  separate  reactants  meet  in  stoichiometric  proportions, 
with  the  temperature  playing  no  role  in  this  balance.  One  may  there¬ 
fore  attribute  the  slight  variation  in  flame  position  between  the  con¬ 
stant  and  variable  cases  to  the  dependence  of  species  distribution  on 
the  overall  density.  As  with  the  flame  temperature,  the  predicted  val¬ 
ues  of  Xf  at  extinction  by  the  two  formulations  is  very  close,  with  the 
difference  being  less  than  0.1%. 

In  Fig.  10,  we  show  the  degree  of  reactant  leakage  through  the 
reaction  zone  for  the  constant  and  variable  density  formulations.  In 
both  cases  the  mass  fractions  tend  to  zero  at  low  strain  rates,  as  it 
should.  Elsewhere,  the  variable  density  solution  shows  a  slight  devi¬ 
ation  in  the  leakage  of  Y3  and  Y2  from  the  constant  density  solution. 
This  is  more  evident  for  titanium,  because  the  extent  of  titanium  dif¬ 
fusion  is  much  smaller  than  that  of  boron  causing  a  sharper  response 
to  density  variations.  The  difference  between  the  two  solutions  at  ex¬ 
tinction  is  nearly  0.005  for  Y2  and  0.05  for  Y] ,  approximately  10%  of 
the  actual  extinction  values. 

5.  Conclusions 

A  multi-component  mixture  theory  is  used  to  describe  condensed 
phase  diffusive  combustion,  in  particular  for  titanium-boron  reac¬ 
tion  in  a  counterflow  geometry.  The  Maxwell-Stefan  model  for  multi- 
component  mixtures  is  informed  by  the  binary  diffusivities  of  any 
pair  of  species  comprising  the  mixture.  Because  of  the  complications 
arising  from  the  direct  use  of  the  Maxwell-Stefan  model  in  combus¬ 
tion  studies,  a  generalized  Fick  law  is  often  adopted  with  coefficients, 
referred  to  as  Fick  diffusivities,  which  are  often  assumed  to  be  con¬ 
stant.  These  coefficients,  however,  are  concentration-dependent  and 
could  be  related  to  the  binary  diffusivities  via  cumbersome  expres¬ 
sions.  A  significant  contribution  of  this  work  is  the  derivation  of  a 
multi-component  diffusion  model  for  a  three-component  mixture, 
which  reduces  to  the  gaseous  model  when  molecular  weights  and 
diffusivities  for  each  species  are  equal. 

The  proposed  diffusion  expressions  can  be  simplified  by  assum¬ 
ing  equal  molecular  weights.  A  comparison  between  the  two  diffu¬ 
sion  models,  with  equal  and  unequal  molecular  weights,  has  been 
carried  out  in  the  constant  density  limit  (for  simplicity  of  computa¬ 
tions).  The  response  curves  constructed  numerically  show  that  the 
two  diffusion  models  lead  to  identical  results  at  large  strain  rates  as 
expected,  and  vary  by  less  than  2%  near  extinction.  This  suggests  that 
one  could  safely  use  the  constant  molecular  weight  approximation  in 
the  proposed  diffusion  model  with  minimal  loss  in  accuracy,  while 
saving  computational  time  and  complexity. 

Low  temperature  measurements  of  the  boron-titanium  diboride 
binary  diffusivity,  ®23,  fall  in  the  range  of  1CU13  to  10~20  m2/s  while 
the  titanium-titanium  diboride  diffusivity,  f^13,  is  three  orders  of 
magnitude  lower.  These  diffusivities,  however,  are  unavailable  at  the 
combustion  temperature  of  Ti-B  nano-composites,  measured  near 
3000  K.  The  derived  expressions  of  flame  position  and  tempera¬ 
ture  from  the  fast-chemistry  asymptotic  analysis  have  been  used 
to  estimate  the  aforementioned  diffusivities  in  a  manner  which  is 
consistent  with  macroscopically  observed  adiabatic  flame  tempera¬ 
tures.  These  estimates  may  hold  in  more  general  circumstances  be¬ 
cause,  as  noted  from  our  simulations  the  flame  temperature  drops 
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only  slightly  over  the  entire  range  of  strain  rates  sustaining  steady 
combustion. 

Combustion  characteristics  at  the  interface  of  two  condensed 
phase  reactant  streams  that  form  an  opposed  counterflow  are  rep¬ 
resented  by  a  typical  S-shaped  response  curve.  In  this  paper  we  only 
address  the  portion  of  the  S-curve  that  extends  from  a  vigorous  burn¬ 
ing  state,  where  the  flame  temperature  reaches  its  maximum  value 
(the  adiabatic  flame  temperature)  and  both  titanium  and  boron  are 
completely  consumed,  down  to  an  extinction  state  associated  with 
incomplete  combustion  and  with  a  flame  temperature  significantly 
below  the  adiabatic  value.  The  extinction  state  is  identified  when  the 
strain  rate  reaches  a  critical  value,  and  for  larger  strain  rates  a  bal¬ 
ance  between  advection,  diffusion  and  reaction  is  not  possible  un¬ 
der  steady  conditions.  The  other  portion  of  the  S-curve,  which  ex¬ 
tends  to  exceedingly  large  values  of  the  strain  rate  is  associated  with 
chemically  frozen,  or  weakly-burning  states.  Our  analysis  therefore 
suggests  that  in  composite  materials  the  distribution  of  local  reac¬ 
tion  sites  and  local  frozen  sites  would  depend  on  the  local  strain  rate 
distribution. 
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The  relation  between  the  Fick  diffusivities  By  and  binary  diffusiv- 
ities  %j  for  a  ternary  mixture  are  given  by 
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Condensed  phase  energetic  materials  include  propellants  and  explosives.  Their  detonation 
or  burning  products  generate  dense,  high  pressure  states  that  are  often  adjacent  to 
regions  that  are  at  vacuum  or  near-vacuum  conditions.  An  important  chemical  diagnostic 
experiment  is  the  time  of  flight  mass  spectroscopy  experiment  that  initiates  an  energetic 
material  sample  via  an  impact  from  a  flyer  plate,  whose  products  expand  into  a  vacuum. 
The  rapid  expansion  quenches  the  reaction  in  the  products  so  that  the  products  can  be 
differentiated  by  molecular  weight  detection  as  they  stream  past  a  detector.  Analysis  of 
this  experiment  requires  a  gas  dynamic  simulation  of  the  products  of  a  reacting  multi- 
component  gas  that  flows  into  a  vacuum  region.  Extreme  computational  difficulties  can 
arise  if  flow  near  the  vacuum  interface  is  not  carefully  and  accurately  computed.  We 
modify  an  algorithm  proposed  by  Munz  [1],  that  computed  the  fluxes  appropriate  to  a 
gas-vacuum  interface  for  an  inert  ideal  gas,  and  extend  it  to  a  multi-component  mixture 
of  reacting  chemical  components  reactions  with  general,  non-ideal  equations  of  state.  We 
illustrate  how  to  incorporate  that  extension  in  the  context  of  a  complete  set  of  algorithms 
for  a  general,  cell-based  flow  solver.  A  key  step  is  to  use  the  local  exact  solution  for  an 
isentropic  expansion  fan,  for  the  mixture  that  connects  the  computed  flow  states  to  the 
vacuum.  Regularity  conditions  (i.e.  the  Liu-Smoller  conditions)  are  necessary  conditions 
that  must  be  imposed  on  the  equation  of  state  of  the  multicomponent  fluid  in  the  limit  of 
a  vacuum  state.  We  show  that  the  Jones,  Wilkins,  Lee  (JWL)  equation  of  state  meets  these 
requirements. 

©  2015  Elsevier  Inc.  All  rights  reserved. 


1.  Introduction 

Condensed  phase  energetic  materials  include  propellants  and  explosives.  They  are  usually  composed  of  a  mixture  of 
granular  solids  that  include  explosive  or  oxidizing  crystallites,  various  metal  powders  like  aluminum,  sometimes  carbon 
black,  resins  and  plastics.  The  performance  of  the  aggregate  composite  depends  on  the  chemistry  and  the  mechanisms 
of  energy  release,  which  occur  in  nearly  all  phases  of  materials,  gas,  liquid  and  solid.  Propellant  and  explosives  reactive 
decomposition  produces  huge  volume  expansion.  The  products  start  at  near  solid  densities  and  at  high  pressures  and  expand 
to  very  low  densities  and  lower  pressures.  In  the  case  of  explosives,  the  pressure  drops  from  hundreds  of  kilo-bars  to  1 
atmosphere  or  less,  which  is  5  to  6  orders  of  magnitude  across  a  reaction  zone  that  is  often  no  more  than  1/10  of  a 
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Fig.  1.  The  schematic  diagram  of  numerical  simulation  of  initiation,  detonation  and  expansion  for  an  explosive  sample  in  the  TOFMS  experiment. 


millimeter  thick.  This  enormous  pressure  gradient  provides  the  means  to  cut  materials  or  drive  surrounding  materials  to 
large  velocities  by  virtue  of  this  expansion  power.  Likewise,  solid  propellants  in  rocket  motors  vent  gases  with  pressures 
that  range  from  f  to  200  atmospheres,  down  to  vacuum  conditions.  Many  orders  of  magnitudes  of  pressure  change  are 
realized  and  the  pressure  gradient  in  propellant  exhaust  stream  provides  the  means  to  generate  thrust.  The  problem  of 
computing  the  transition  of  material  states  from  very  high  pressure  to  vacuum  or  near  vacuum  states  is  a  generic  one  for 
any  multi-material  simulation  where  two  materials  may  collide  and  have  individual  or  shared  boundaries,  for  which  one 
of  the  materials  is  adjacent  to  a  region  with  very  low  pressure  and  density  or  a  vacuum.  Difficulties  generally  arise  if  the 
region  of  expansion  between  the  high  and  low  pressure  regions  is  not  accurately  computed. 

A  time  of  flight  mass  spectrometry  (TOFMS)  experiment  by  Fossum  et  al.  2],  initiated  small  samples  of  energetic  ma¬ 
terials  by  laser  flyer  plate  impact.  After  impact,  the  products  expanded  into  a  vacuum  region  where  spectroscopic  analysis 
of  the  products  was  performed.  A  larger  version  of  a  similar  experiment  was  carried  out  by  Blais  et  al.  [3]  that  used  a 
large  quantity  of  explosive  to  drive  an  ampule  of  nitromethane,  which  detonated  upon  being  shocked.  The  premise  of  the 
TOFMS  experiments  is  that  the  pressure  drop  caused  by  the  expansion  in  a  long  vacuum  region  freezes/quenches  the  reac¬ 
tion  amongst  product  species  in  a  sequential  manner,  ordered  by  the  events  in  the  reaction  zone  in  the  material  that  was 
established  prior  to  the  expansion.  In  order  to  interpret  TOFMS  experiments,  accurate  model  simulations  of  the  gas  dynam¬ 
ics  of  reacting  multi-component  mixtures  are  required.  The  range  of  thermodynamic  states  is  extreme.  The  reactants  start 
at  standard  room  pressure  and  temperature  conditions,  are  raised  to  high  pressure,  density  detonation  or  highly  shocked 
states  in  condensed  materials,  convert  to  (mostly)  gas  species  products  that  expand  to  a  vacuum  or  near-vacuum  state.  Such 
simulations  are  very  challenging  to  carry  out  since  non-ideal  equations  of  state  forms  must  be  used  for  the  reactants  and 
products,  and  the  location  of  the  vacuum/materials  boundary  must  be  calculated  in  a  precise  manner. 

Computational  methods  and  techniques  are  not  widely  available  for  multi-component  reactive  flow,  when  the  compo¬ 
nents  are  subjected  such  huge  ranges  in  the  thermodynamics  states.  Since  we  had  an  interest  in  finding  an  accurate  and 
robust  way  to  simulate  the  TOFMS  experiments,  we  decided  to  make  a  modification  to  the  basic  vacuum  interface  tracking 
algorithm  pioneered  by  Munz  [1],  Our  extension  is  used  in  combination  with  a  (now)  standard,  higher  order,  cell-based, 
totally  variation  diminishing  (TVD)  Euler  scheme,  in  a  fairly  general  multicomponent  framework.  One  should  be  able  to 
use  our  method  to  compute  the  approach  to  the  vacuum  in  the  continuum  limit,  and  combine  it  with  simulations  of  the 
Boltzmann  equation  or  with  molecular  dynamic  simulations.  The  continuum  simulations  we  describe  would  generate  the 
near  vacuum  continuum  flows  as  a  far-field  limit,  or  be  used  to  establish  averaged  initial  conditions  for  molecular  based 
simulations.  In  this  paper,  we  describe  these  algorithms  and  implementations  in  detail,  and  present  worked  examples  that 
simulate  the  flow  in  the  TOFMS  experiment  as  a  targeted  application. 

Fig.  1  shows  a  sketch  of  a  typical  sequence  of  events  that  occur  in  a  simulation  of  the  quenching  of  blast  products  that 
expand  into  a  vacuum  region  on  the  right,  a)  The  explosive  sample  and  vacuum  regions  are  initially  separated,  prior  to 
impact  by  the  flyer,  b)  If  the  sample  detonates,  a  detonation  shock  and  its  supporting  reaction  zone  propagate  through  the 
explosive,  c)  The  detonation  wave  hits  the  vacuum  interface,  d)  Then  a  rarefaction  travels  back  into  the  detonation  reaction 
zone  structure.  As  the  pressure  and  density  drop,  the  reactions  amongst  product  species  slow  as  the  products  flow  into  the 
vacuum  section,  past  the  detector  which  monitors  the  mass  concentration  of  the  products  at  a  fixed  probe  locations.  Since 
the  number  density  drops  in  the  vacuum  region,  the  species  are  thought  not  to  undergo  significant  collisions  and  hence  no 
further  chemical  changes  occur. 
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1.1.  Summary  of  the  algorithms 

The  original  work  of  Munz  [1]  proposed  an  approximation  method  that  computes  the  fluxes,  appropriate  to  the  gas- 
vacuum  interface,  and  tracked  the  motion  of  the  gas-vacuum  boundary  for  an  inert  ideal  gas.  Our  extended  algorithm 
allows  for  reactive  flow  with  general,  non-ideal  equation  of  state  (EOS)  in  Mie-Gruneisen  (MG)  form  for  a  reacting  multi- 
component  mixture,  adjacent  to  a  vacuum  interface.  The  formulation  can  be  used  to  include  an  arbitrary  number  of  chemical 
components  or  other  state  variables  like  porosity  or  particulate  density,  that  could  be  included  in  the  constitutive  descrip¬ 
tion.  Our  standard  multi-component  flow  model  has  only  one  material  velocity  for  the  mixture.  The  corresponding  gas 
dynamics  formulations  have  only  forward,  backward  acoustic  characteristic  and  one  material  characteristic.  The  exact  form 
of  the  vacuum  Riemann  problem  can  thus  be  handled  in  a  standard  and  familiar  way,  albeit  with  complexities  associated 
with  the  EOS  that  we  describe  in  some  detail.  We  note  that  Lee  et  al.  [4]  recently  considered  the  Riemann  problem  for  the 
MG  EOS  in  their  analysis  of  difficulties  associated  with  contact  discontinuities. 

The  extended  vacuum  tracking  algorithm  is  described  in  the  context  of  other  supporting  methods  and  solver  algorithms. 
We  use  a  cell-based  reactive  Euler  solver  that  employs  local  cell-based  Riemann  problems  for  general  EOS  to  compute  the 
cell  edge  fluxes.  Roe-average  fluxes  are  used  in  the  manner  originally  developed  by  Glaister  [5],  and  extended  by  Xu  [6],  and 
Xu,  Aslam  and  Stewart  [7]  for  chemically  reacting  species.  The  reactive  Euler  equations  are  solved  by  source  splitting  with 
the  first  update  as  a  solution  to  the  Euler  equations  without  source  term,  followed  by  an  update  of  the  (chemical)  source 
terms.  We  will  reference  to  the  first  update  as  the  hydrodynamic  update  and  the  second  update  as  the  source  term  update.  We 
have  used  5th  order  Weighted  Essentially  Non  Oscillatory  (WENO)  with  3rd  order  Runge-Kutta  (RI<)  for  the  hydrodynamic 
update.  Updates  to  the  location  and  velocity  of  vacuum  interface  of  material,  and  the  special  cell  boundary  flux  calculations, 
use  the  extensions  to  the  Munz  [1]  tracking  scheme.  We  have  verified  our  schemes  against  various  examples  found  in  the 
literature  [1,5,8], 

Care  must  be  taken  to  accurately  integrate  the  chemical  source  terms,  because  the  rates  in  most  kinetics  schemes  widely 
vary  and  can  have  both  fast  and  slow  time  scales  that  are  disparate.  Colella  et  al.  [9  were  amongst  the  first  to  point  this  out 
for  Arrhenius  kinetics,  applied  to  the  calculation  of  gas-phase  reacting  flow.  In  this  case  it  was  advantageous  to  analytically 
integrate  an  approximate  form  of  the  source  terms  across  a  time  step.  However  source  term  splitting  treatments  have  been 
proven  to  work  well  (dating  back  to  Colella’s  original  work)  and  are  known  to  have  good  stability  properties.  In  his  thesis 
Xu  [6]  demonstrated  that  an  algorithm  of  the  same  type  we  use  here  for  non-ideal  EOS  forms,  was  second  order  accurate. 

If  one  uses  an  inaccurate  numerical  method  for  the  basic  hydrodynamic  solver,  significant  problems  arise  from  inaccu¬ 
racies  in  the  evaluation  of  the  sound  speed  c  =  ^/dp/dp\s  in  the  vicinity  of  the  vacuum  region,  in  which  case  tracking  the 
vacuum  material  interface  becomes  difficult  or  breaks  down.  Therefore  in  order  to  compute  the  flux  state  and  velocity  of 
the  vacuum  material  interface,  we  work  out  the  exact  form  of  the  isentropic  expansion  fan  that  vents  to  the  interface.  This 
connects  the  vacuum  state  to  the  state  in  the  adjacent  cell,  and  thus  supplies  the  combined  method  with  a  physically  accu¬ 
rate  value  for  the  vacuum  interface  velocity.  Since  we  allow  for  general  equation  of  state,  we  must  insist  that  the  limiting 
behaviors  of  the  EOS  forms  are  such  that  the  vacuum/material  interface  problem  is  well-posed.  Fortunately,  the  conditions 
that  define  a  well-posed  vacuum  Riemann  problem  (VRP)  have  been  considered  by  Liu  and  Smoller  [10]  and  they  provide 
criteria  on  limiting  functional  form  of  the  EOS  as  the  material  expands  to  the  vacuum  state  of  zero  mass  density.  These 
conditions  must  be  incorporated  into  the  vacuum  interface  tracking  method.  We  show  that  the  Jones,  Lee,  Wilkins  (JWL) 
EOS  form,  commonly  used  for  product  components,  satisfies  these  conditions. 

1.2.  Paper  outline 

Section  2  describes  basic  model  formulations  that  include  the  conservation  laws,  the  mixture  equation  of  state  for  a 
multi-component  mixture  with  commonly  used  EOS  forms  for  the  components,  that  include  ideal  gas  EOS  forms  and 
Mie-Gruneisen  EOS  forms.  Details  are  given  for  the  JWL  EOS  form  for  reactant  and  product  components.  Details  for  two 
component  mixture  closures  are  given.  One  is  the  standard  pressure-temperature  equilibrium  closure  (PT-EQB),  and  the 
other  is  an  often-used,  simpler,  pressure  equilibrium  closure  (P-EQB),  that  specifies  a  volume  ratio  between  the  components 
(instead  of  enforcing  strict  temperature  equilibrium).  Section  3  describes  the  numerical  methods,  with  a  brief  review  of 
basic  methods  for  the  compressible  Euler  equations.  A  brief  description  of  the  Riemann  solver  for  general  EOS  is  given  fol¬ 
lowed  by  the  detailed  discussion  of  the  implementation  of  the  vacuum  tracking  method  for  a  general  EOS  in  a  reactive  flow. 
Section  4  presents  examples  and  applications.  The  first  example  describes  a  condensed  explosive  with  ideal  gas  EOS  forms 
for  both  the  reactants  and  products.  This  is  followed  by  a  detailed  simulation  of  a  facsimile  of  the  TOFMS  experiment  where 
a  thin  sample  of  PETN  is  hit  by  an  impactor  and  whose  products  are  vented  into  the  long  vacuum  region  of  the  experiment. 
The  components  are  modeled  with  non-ideal  MG  EOS  forms.  We  also  discuss  the  sensitivity  of  the  simulation  results  to  the 
different  mixture  closure  conditions,  pressure,  temperature  equilibrium  (PT-EQB)  and  pressure  equilibrium  (P-EQB). 

2.  Basic  model  formulations 

Here  we  give  the  formulations  for  the  physical  models  of  interest.  The  models  apply  to  a  mixture  of  compressible  fluids 
with  N  distinct  species.  The  fluid  mixture  is  assumed  to  have  a  single  mixture  material  velocity  u.  pressure  p,  and  generally 
a  single  temperature  T  (although  that  is  not  strictly  required).  Each  species  will  be  assumed  to  be  specified  by  its  separate 
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equation  of  state,  and  the  mixture  energy  that  is  used  to  define  the  equilibrium  potentials  will  neglect  mixing  terms,  that 
can  sometime  be  included  in  the  theory  of  liquid  mixtures.  Hence  the  mixture  is  ideal,  although  the  equation  of  state 
(EOS)  forms  are  not  necessarily  those  of  the  ideal  gases.  We  assume  that  the  EOS  forms  provided  by  the  model  allow  for 
a  continuous  and  smooth  approach  to  zero  density,  i.e.  to  the  vacuum  state.  That  will  allow  us  to  calculate  reasonable 
mixture  velocities  and  track  the  location  of  vacuum  material  interface.  It  is  understood  that  the  continuum  limit  breaks 
down  at  a  true  vacuum  and  the  model  only  predicts  some  limiting  approximation  (from  the  continuum  limit  side)  of  the 
states  near  the  material/vacuum  interface.  The  description  of  the  transition  to  vacuum  on  the  molecular  scale  requires  the 
implementation  of  additional  models,  that  are  not  the  subject  of  this  work.  Thus  the  flow  simulations  as  shown  in  the 
section  for  examples  are  restricted  to  extremely  low  but  non-zero  densities,  where  the  continuum  approximations  are  still 
reasonable. 


2.1.  The  governing  equations 


It  is  useful  to  present  all  three  forms  of  the  governing  equations,  the  primitive,  conservation  and  characteristic  forms. 
The  basic  numerical  scheme  in  the  body  of  the  fluid  uses  the  conservative  form,  while  the  vacuum  interface  approximation 
require  the  characteristic  form.  The  primitive  form  of  the  governing  conservation  laws  for  mixture  mass,  momentum,  energy, 
and  mass  fraction  of  identified  chemical  species  are  respectively 


p  +  pV  •  v  =  0, 

(1) 

p  v  =  -Vp , 

(2) 

pe  +  pV  •  v  =  0, 

(3) 

pU  =  rXi , 

(4) 

where  we  use  the  “dot"  notation  is  used  for  the  material  derivative  ()  =  3/3 1  +  v  •  V.  The  variables  p  =  1/v,  e,  v  and 
Aj  are  the  mixture  density  (its  reciprocal,  specific  volume),  specific  internal  energy,  mixture  velocity  and  mass  fraction  of 
the  ith-species  respectively.  The  one-dimensional  conservative  form  of  the  Euler  equations  for  reactive  flow  in  the  multi- 
component  mixture  (illustrated  for  two  independent  components/reactions,  A.i,  X2)  are: 


and 


3U  3F  (U) 
~dt  +  dx 


=  S(U) 


-  p  - 

pu 

-  0  - 

pu 

pu2  +  p 

0 

U  = 

pE 

F  = 

( pE  +  p)u 

S  = 

0 

px  i 

puX  i 

pt\ , 

-px  2- 

_  puX 2  - 

-prX2- 

1 

E  =  e  - 1 — u 
2 


2 


(5) 


(6) 


(7) 


Riemann  problems  and  the  methods  for  solving  these  hyperbolic  equations  are  based  on  the  characteristic  form  of  the 
equations  which  we  list  below  in  one  spatial  dimension,  x 


3 P  ,  ,  .  dp 

-  +  IU  +  C-  +  PC 

dp  dp 

- h  (u  —  C) - pc 

dt  dx  P 


du  du' 

dt  ’  dx 

du  du' 

- b  (U  —  C) - 

dt  dx 


=  pc\a  r), 


=  pc2  (a  ■  r) , 


Tr+u^Z-c2  1=  Pc^°  ■  r) 


dt 


dx 


dt 


dx 


dki  dXi 
~df  +  u~dx  ~ rXi ' 


The  sound  speed  c2  and  the  thermicity  a  ■  r  are  given  by 


c2=  (§1  P*+P)VZ 


(  ff  I  v,A.,  ) 


with  Oi  = 


1  (de/dki)\p,v 
pc 2  (de/dp)\v,Xi 


and  a  ■  r  =  ^  airXi  , 


(8) 

(9) 

(10) 

(11) 


(12) 


respectively,  and  where  the  specific  internal  energy  is  of  the  form  e(p,  v,A.,-).  The  equations  above,  written  as  ordinary 
differential  equations  on  their  respective  characteristics,  are  given  by 
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dp 

du 

9 

dx 

(13) 

dt 

+  fc2i 

=  pc1  (a  ■  r ) 

on  C+: 

—  =  u  +  c, 
dt 

dp 

du 

dx 

(14) 

dt 

~PCdt 

=  pc\a  ■  r) 

on  C_: 

—  =  u  —  c , 
dt 

dp 

dt 

dt 

=  pc2  (a  ■  r) 

on  C0: 

dx 

dt=U ’ 

(15) 

dXi 

dx 

(16) 

dt 

-r»  onc°:  *= 

u . 

2.2.  Mixture  equation  of  state 


We  assume  that  the  fluid  mixture  is  comprised  of  N  chemical  components.  A  complete  equation  of  state  is  supposed  to 
be  given  for  each  component  (species).  One  can  suppose  that  a  Gibbs  free  energy  form  is  provided  for  each  component  that 
specifies  that  energy  at  a  fixed  temperature  and  pressure  in  the  form  gi(p,  T ),  which  is  a  complete  equation  of  state,  Callen 
[11],  from  which  we  can  write  the  specific  internal  energy  and  pressure,  volume  temperature  equation  of  states,  in  the  form 

e,(p,Vi)  and  v,- =  v,-(p,  T) .  (17) 

The  internal  energy  and  volume  are  extensive  thermodynamic  quantities,  so  the  contribution  from  each  component  add 
in  a  fixed  volume  of  mixture  in  proportion  to  the  relative  masses  of  the  components.  Additional  contributions  due  to 
the  intermolecular  mixing  are  neglected.  If  the  mass  fraction  of  each  component  is  X,,  then  the  specific  energies  and  the 
volumes  of  the  mixture  can  be  written  as  the  sums 


N 

e(p,  T ;  7,)  =  ^e,(p,  v,(p,  T))Xj  and 

i=l 


N 

v(p,T;Xi)  =  '^vi(p,T)Xi, 

i=i 


(18) 


with  Xi  =  1.  For  fixed  mixture  energy  and  volume  and  internal  composition,  the  first  and  second  equation  of  Eq.  (18) 
are  implicit  equations  for  the  pressure  and  temperature,  (p,  T). 

We  require  the  regularity  conditions  on  the  equation  of  state,  called  the  Liu-Smoller  conditions  (see  [10]),  when  the 
density  approaches  to  zero  density  on  an  isentrope.  The  pressure  relations  modeled  by  the  mixture  EOS  (18)  should  satisfy 


P(0)  =  0, 


9p(0) 

dp 


Is  =  0, 


and 


>0, 


92p 

— ^-]s>0  for  p  >  0 . 
dp 2 


(19) 


In  order  to  model  a  particular  energetic  or  explosive  material,  one  must  identify  all  the  components  and  their  EOS  forms 
for  each  component.  The  ideal  gas  EOS  form  is  often  used  in  detonation  theory,  [12],  and  is  also  the  EOS  form  that  was 
used  by  Munz  to  establish  his  original  vacuum  tracking  algorithm.  We  describe  those  common  EOS  forms  next. 


2.2.1.  Ideal  gas  EOS  forms  for  the  components 

When  all  the  components  are  specified  by  the  ideal  EOS  form  with  constant  specific  heats  of  formation,  we  write 


1 

e;  = - -pvj  +ei0, 

Vi  ~  1 


(20) 


where  y,-  =  cpi- /cvi  (the  ratio  of  specific  heats)  and  e;o  is  the  heat  of  formation  of  the  i-th  component.  We  also  use  the 
notation  &),-  =  (y,-  —  1).  The  component  gas  constant  is  Rj  =  cV!(yj  —  1),  and  the  p,  v,-,  T  EOS  is  given  by 


pvt  =  RjT 


(21) 


The  mixture  energy  and  mixture  volume  are  given  by 

e  =  p  V  — Xi  +  V  e,0  Xi ,  and  v  =  -  V  Rt  Xj . 
Vi  —  1  t—1  p  x—1 

i  i  i 

The  first  part  of  Eq.  (22)  can  also  be  written  as 


(22) 


(23) 


2.2.2.  A  Mie-Gruneisen  EOS  form  for  components 

A  e(p,  v)  form  of  Mie-Gruneisen  (MG)  EOS  for  components  is  explicitly  linear  in  the  pressure  p  and  generally  nonlinear 
in  the  specific  volume,  and  can  be  written  (dropping  the  i-subscript  temporarily)  as 

e(p,  v)  =  es(v)  +  ^(p  -  ps(v))  =  A(v)  +  C(v)p 


(24) 
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The  functions  es(v),  ps(v),  and  v/T  are  functions  that  are  determined  from  EOS  calibration  experiments.  Often  es(v)  and 
ps(v )  are  identified  as  the  internal  energy  and  pressure  at  reference  states  on  some  specified,  well-identified  and  measur¬ 
able  thermodynamic  process,  such  as  the  states  on  a  p,  v  expansion  isentrope,  or  a  shock  Hugoniot.  Here  they  are  simply 
regarded  as  known  functions  of  the  specific  volume,  v.  The  Mie-Gruneisen  coefficient  r  is  often  assumed  to  be  a  constant 
or  simply  a  function  of  v.  We  note  that  the  ideal  EOS  forms  are  in  this  MG  class,  so  the  MG  forms  are  inclusive  of  the  ideal 
EOS  forms. 


2.2.3.  The  Jones,  Lee,  Wilkins  (JWL)  EOS  form  for  components 

The  Jones,  Lee,  Wilkins  (JWL)  EOS  form  for  components  is  in  the  MG  class,  and  is  probably  the  most  commonly  used 
EOS  form  for  engineering  studies  of  reactive  flow  in  explosives.  For  example,  the  chemical  equilibrium  code  Cheetah  [13] 
automatically  generates  JWL  equation  of  state  for  reaction  products.  The  JWL  EOS  forms  are  often  paired  with  the  “Ignition 
and  Growth”  reaction  rate  models  that  were  pioneered  by  Lee  and  Tarver  [14].  Many  highly  used  explosive  compounds  such 
as  HMX,  RDX,  PETN,  TATB,  TNT,  etc.  are  modeled  with  MG  EOS  forms,  paired  with  Ignition  and  Growth  kinetic  rate  models 
that  describe  the  rate  of  overall  decomposition  from  reactants  to  products. 

A  standard  convention  is  to  use  a  lower  case  letter  for  specific  internal  energy  (as  in  e)  and  a  upper  case  letter  for  energy 
per  unit  volume  E.  An  upper  case  V  to  represent  the  specific  volume  ratio  relative  to  a  reference  volume  (as  in  V  =  v/vo). 
Then  for  each  component  the  JWL  and  thermal  EOS  is  given  by 


p  =  At  exp-*1'1''  +B,  exp  R2iv>  +^T 


with 


Vi 

Vt  =  —  , 
Vo 


(25) 


and  the  JWL  EOS  form  are  given  by  the  e(p,  v,-)  mechanical  EOS  given  by 


Ei  =  £f(V{)  +  —  [p  -  P- (V,)]  +  E0i , 

COi  L  J 


(26) 


with 


Ej(v)  = 


Ae 


-Rlivi  JJe~R2ivi 
- -  + 


R2i 


j  and  psi(Vi)  =  (^Aie~RliVi +  Bie~RaV‘Y 


(27) 


With  e,  =  E,  vo  one  has 

ei  =  e- (vO  +  —  [p  -  p- (Vj)]  +  e0i- ,  (28) 

a>i 

with  e?(v,)  =  E^vq.  The  constants  A,,  B,-t  Ru,  R 2i,  COi  and  Cvi  =  cVi/vo  are  the  JWL  component  EOS  constants;  eoi  is  the 
reference  energy  and  vq  is  a  standard  reference  volume  of  the  original,  unshocked  mixture. 


2.2.4.  A  mixture  EOS  of  reactant  and  product  components  for  the  JWL  EOS  forms 

Next  suppose  we  specify  that  the  mixture  is  composed  of  only  two  species/components  reactants  (R)  and  products  (P). 
Both  are  assumed  to  be  modeled  by  the  JWL  EOS  forms.  The  material  decomposes  by  the  reaction  R  —>■  P,  during  which  the 
reactant  and  product  are  assumed  to  be  in  pressure  and  temperature  equilibrium.  The  product  mass  fraction  is  taken  be  X 
while  the  reactant  mass  fraction  be  1  —  X.  An  R  subscript  is  used  to  denote  reactant  and  P  for  product.  Given  the  energy  e 
and  a  mixture  volume  v,  the  constraints  (18)  become 

e  =  XeP(p,  vP(p,  T))  +  (1  -X)eR(p,  vR(p,  T)),  (29) 

v=Xvp(p,T)  +  (L  -X)vR(p,T).  (30) 

Since  the  JWL  forms  (25)  and  (26)  are  not  generally  linear  in  the  specific  volume,  the  volume  dependence  is  implicit 
in  p,  T.  Given  specified  values  for  e,  v,  X,  Eqs.  (30)  and  (29)  are  two  equations  for  p  and  T.  The  solution  completes  the 
EOS  description  for  the  mixture,  and  defines  the  PT-EQB  closure.  There  are  many  ways  to  find  the  roots  but  all  involve 
an  iteration  to  find  roots  of  a  pairs  (or  n-tuples)  of  equations.  For  the  specific  choice  of  the  MG,  JWL  EOS  forms  for  two 
components,  the  problem  of  finding  the  (p,  T)  root  pairs  can  be  reduced  to  a  single  scalar  unknown.  Next  we  present  a 
simple  robust  algorithm  that  can  be  used  to  find  the  roots. 

2.2.5.  Determination  of  the  ( p ,  T)  (Gibbs)  equilibrium  states 

We  present  a  simple  way  to  compute  the  PT-EQB  for  JWL  EOS  forms,  that  can  be  used  for  N  components;  illustrated  for 
three  that  can  be  extended  simply  to  N.  The  mechanical  and  thermal  EOS  forms  can  be  represented  as 

e,(p,Vi)  and  p  =  p(v;,T), 

for  some  fixed  temperature  and  pressure,  where  e,-,  v,  are  the  specific  energy  and  volume  for  that  component.  The  mixture 
energy  and  volume  are  represented  as  a  sum  over  the  components 

e  =  X\et  +  >1262  +  ^3^3  ,  and  v  =  Xj  Vi  +  X2V2  +  ^3 V3  . 


(31) 


164 


S.  Choi  et  al.  /  Journal  of  Computational  Physics  296  (2015)  158-183 


Use  Vi  as  a  reference  volume  and  define  the  component  volume  ratios 
<t>2  =  V2/V1  and  <t>3  =  v3/vi, 


(32) 


which  combined  with  the  second  equation  in  Eq.  (31)  solves  for  the  component  volumes  in  term  of  the  mixture  volume 
and  the  ratios,  to  obtain 


Vi  = 


V2  = 


vt>2 


+  X2^2  "F"  ^-3<I,3  A}  +  X2&2  “b  7.3<t>3 

Using  the  forms  (28)  for  e,-  in  (31)  obtains 

e  =  p  V  —  A.,-  +  V'[e'(Vi)  -  —  P- (Vi)]Xj  +  Ve.oA;  ■ 

COj  >  &>■  2__^ 

2  2  2 


V3  = 


vd>3 


X\  +  X2Q2  “h  ^<3^3 


(33) 


(34) 


Solving  for  p  gives 


P  = 


e  -  y'[e?(vi)  -  —  p- (vOJXi  -  V  ei0Xi 

u>i 


/iE 


COi 


(35) 


The  pressure  is  now  explicitly  determined  by  direct  evaluation  once  the  mixture  energy  e  and  the  volume  ratios  $2.  $3 
are  known.  Thus  one  needs  two  (2)  residuals  to  determine  the  correct  values  of  <l>2  and  $3.  Those  residuals  come  from 
setting  the  temperature  of  the  component  (species)  equal.  For  the  case  of  three  components  there  are  two  such  equivalences, 
i.e.  T  =  T 1  =  T 2,  and  T 2  =  T3.  For  N  components  there  are  N  — 1  volume  ratios  and,  likewise  N  —  1  temperature  equivalences. 
For  the  3-component  case,  with  the  JWL,  MG  forms,  these  equivalences  are 


T 

T 


—  1 

1=1^1 

\p-(A2e-^V2  +  B2e-R22vA] 

CO  1  1 

L  V  /  J 

Cl  a> 2  1 

L  V  /  J 

V3 1 

\p-(A3e-Rlvi  +  B3e-RlvA] 

1  =  1^1 

\p-(A2e-RiV2+B2e-R22vA] 

(D\  \ 

L  V  /  J 

Cl  co 2  1 

L  V  /  J 

recast  as  the  pair 


J_V± 

Cl  coi 


[p-Pi(Vi)] 


J_V2 
cl  U>2 


[p-ps2(v2)\ 


1  V3 

cl  a>3 


[P-PKV3)] 


_1_  V2 
Cl  U>2 


[P-PS2(V2)] 


(36) 

(37) 

(38) 


that  defines  the  two  residuals.  For  a  two  component  case  only  one  residual  is  required. 

The  roots  are  determined  by  the  following  steps.  1)  Specify  the  mixture  values  of  v  and  e  and  the  reaction  progress 
variables,  X ;.  2)  Make  an  initial  guess  for  $2  and  $3  (say)  and  then  the  values  of  the  component  volumes  (for  the  example: 
Vi,  v2  and  V3  are  computed.  3)  Then  the  residuals  defined  by  Eq.  (38)  can  be  evaluated.  4)  Iterate  using  the  residuals  to  find 
the  root.  Any  reasonable  root  procedure  can  be  used.  Newton-Raphson  works  well  and  converges  in  a  few  steps.  In  special 
cases  direct  searches  can  be  used.  The  initial  guess  for  the  volume  ratios  are  usually  determined  from  the  previous  (time 
step)  physical  state,  so  good  quality  initial  guesses  are  in  abundance.  The  convergence  of  the  root  finding  procedure  then 
determines  all  the  component  ratios,  which  in  turn  can  be  used  to  directly  compute  the  pressure  from  (35)  and  temperature 
from  (36)  or  (37). 


3.  Numerical  methods 


3.1.  The  Euler  solver 


A  typical  numerical  solution  of  our  equations  would  use  operator  splitting  that  updates  a  discrete  approximation  to  the 
source  free  Euler  equations,  followed  by  an  update  of  the  source  terms.  The  latter  generally  requires  a  stiff  ODE  solver  that 
can  handle  the  stiff  reaction  terms.  In  our  implementation  we  used  the  splitting  methods  described  by  Toro  [8].  A  Roe-type 
solver  is  used  to  update  the  source  free  Euler  equations,  and  a  stiff  ODE  package,  Livermore  Solver  of  Ordinary  Differential 
Equations  (LSODE)  [15],  is  used  to  update  the  source  terms  for  interior  points  in  the  domain.  The  first  update  takes  the 
initial  data  at  the  current  time  tn  and  solves  the  source  free  Euler  equations  and  returns  the  updated  states  at  tn+1 .  Thus 

CYCLE  1:  The  Hydrodynamic  update:  Solve 

Ut  +  F(U)X  =  0 ,  (39) 

subject  to  the  initial  condition  U(x,  tn)  =  Un,  and  return  the  solution  at  tn+1  so  that  uhydro  =  U(x,  tn+1). 

The  second  source  term  update  uses  the  first  update  as  the  initial  conditions  at  tn,  and  then  solves  just  the  temporal  ODES 
with  the  source  terms,  to  generated  the  final  update  at  fn+1.  Thus 
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CYCLE  2:  Source  term  update:  Solve 

dU 

^=S(U),  (40) 

subject  to  the  initial  condition  U(x,  tn)  —  U,,ydro,  and  return  the  solution  at  tn+1,  U(x,  tn+1)  =  Un+1. 


The  hydrodynamics  updates  Udydro  obey  (39)  and  their  cell  state  averages  U"+1/,2,  centered  at  i  +  1/2,  are  advanced  in 
the  standard  way  as 

Ui+l1/2=Ui+,/2'  ^[pi+l  -F,l  ,  (41) 

where  F,+1  and  F,-  are  numerical  fluxes  at  the  cell  boundaries.  The  numerical  flux  is  computed  by  using  Roe’s  approximation 
as  follows: 

F;+1  =  1  (F(UL)  +  F(Ur))  H  •  (42) 


The  vectors  and  Ur  in  Eq.  (42)  are  the  state  vectors  at  cell  interface  at  ith  grid  point  and  these  two  vectors  were  obtained 
by  using  WENO  scheme,  which  were  first  introduced  by  Harten  et  al.  (see  [16]).  The  details  of  WENO  reconstruction  are 
found  in  [7,17,18]  for  example.  The  other  quantities  dj,  A,-,  RJ  in  Eq.  (42)  are  the  vector  a,  matrices  functions  A,  and  R 
applied  at  Roe  averaged  state  U.  Those  computations  require  the  Jacobian  of  the  exact  system  given  by 


0 

1 

0 

0 

0 
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0 
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^2 

0 

0 

U 

(43) 


where 


H  =  e+^u2  +  -  and  c2  =  p  p  +  f- p  e ,  (44) 

2  P  P2 

and  the  comma  notation  subscript  denotes  the  partial  derivative  with  respect  to  thermodynamics  argument. 

The  Jacobian  can  be  diagonalized  and  has  real  eigenvalues  (the  characteristic  wave  speeds)  and  real  eigenvectors  accord¬ 
ing  to 

J  =  R  AR1,  (45) 


where  R  is  identified  as  the  matrix  composed  the  right  eigenvectors  of  J  (columns) 


0 
0 

Ea 2 

P,e 

o 

l 

and  A  is  the  diagonal  matrix  corresponding  to  the  eigenvectors 
A  =  diag[(u  —  c),  u,  ( u  +  c),  u ,  u]  . 

The  Roe  averaged  states  are  defined  by 


R(U)  = 


1 

u  —  c 
H  —  uc 

_  ^2 


H-p/p-p 

At 

^2 


1 

u  +  c 
H  +  uc 

^2 


0 

0 

PAi 

P,e 

1 

0 


0=  [p,u,e,X  i,X2], 

where 

P  =  V Pr  Pl, 

~  _  ^/PlUl  +  ->/PrUr 
Vp~l  +  «J~pr 

-/pLeL  +  *J7>ReR 
sfPl  +  VPR 
_  */PL^\L  +  sfPR  ^1R 
■Jpl  +  -J~Pr 

-  sTPL  ^2L  +  VPR  X2r 

A  2  =  - ; - 

■/PL  +  «fPR 


(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 

(53) 
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Head  of  rarefaction  Vacuum  interface 


Fig.  2.  (a)  A  solution  of  Vacuum  Riemann  Problem:  This  solution  is  not  possible  if  UR  is  a  vacuum  state,  (b)  A  possible  solution  of  a  Vacuum  Riemann 
Problem. 


The  enthalpy  H  in  J(U)  is  also  computed  in  the  same  way, 

~  =  i/PL  Hi  +  ^/Pr  Hr  (54) 

VpT + SPr 

Detailed  derivations  of  Roe  average  and  the  pressure  derivatives  of  the  averages  can  be  found  in  Roe  [19,20],  Glaister  [5], 
Xu  [6j  and  elsewhere. 

3.2.  Vacuum  Riemann  problem  (VRP)  and  vacuum  tracking  method  for  non-ideal  EOS 


The  vacuum  Riemann  problem  (VRP)  is  by  definition  a  Riemann  problem  with  vacuum  initial  condition  on  one  side 
and  a  constant  state  at  another  side  of  a  computational  domain.  The  vacuum  is  defined  such  that  the  mass  density  in  the 
geometrical  domain  of  the  vacuum  is  zero,  i.e.  there  is  an  absence  of  material.  Consider  the  situation  with  the  presence  of 
vacuum  in  the  right  state  (x  >  0).  In  the  vacuum  region,  the  conserved  variables  are  zero  so  that  the  sound  velocity  and 
pressure  will  also  be  zero.  The  VRP  for  the  source  free  Euler  equations  can  be  written  as 


9U  9F(U) 

- 1 - -  =0, 

dt  9x 

and  the  initial  condition  is 


(55) 


U(x,  0)  = 


(Pl,  ul,  eL,  Xu,  XL2) 

(0,  0,  0,  0,  0) 


for  x  <  0 , 
for  x  >  0 . 


(56) 


The  vacuum  Riemann  problem  has  a  solution  structure  that  is  different  from  a  conventional  Riemann  problem,  but  it  can 
be  obtained  as  a  limit  of  the  conventional  Riemann  problem  [21,8], 

We  start  with  a  brief  description  of  the  vacuum  Riemann  problem  and  Munz’s  tracking  method  and  then  extend  the 
algorithm  for  a  general  EOS  of  the  form  p  =  p(v,e,X  1,^2)-  Fig-  2(a)  shows  a  Riemann  problem  that  consists  of  constant 
left  state  Ul,  that  is  connected  to  a  left-going  fan,  that  terminates  on  a  constant  state  Ui  that  has  a  particle  velocity  that 
defines  a  contact  line  in  x,  t  space.  The  right  state  Ur  is  connected  by  a  right  going  shock  to  another  constant  state  U2 
whose  pressure  and  particle  velocity  (but  not  necessarily  the  density)  must  match  across  the  contact  discontinuity  to  the 
constant  state  U^  Fig.  2(b),  shows  a  left  state  Ul,  that  is  connected  by  a  single  (left-going)  rarefaction  fan  that  terminates 
on  the  vacuum  state  denoted,  Uo.  In  what  follows,  the  velocity  of  the  material  vacuum  interface  (MVI)  is  uv. 

Next  we  recount  the  argument  that  the  scenario  shown  in  Fig.  2(a),  cannot  be  used  as  a  solution  to  the  VRP  (or  rather 
its  analysis  collapses  the  configuration  to  that  shown  in  Fig.  2(b)).  Given  the  configuration  shown  in  Fig.  2(a),  the  jumps  in 
the  state  across  the  right-going  shock  with  speed  S  (say)  must  obey  the  standard  Rankine-Hugoniot  conditions  for  fixed 
composition, 

P2U2-  PoUo  =  S(p2-  Po),  (57) 

p2uj  +  P2-  {pou20  +  po)  =  S(p2u2  -  pou0) ,  (58) 

U2(P2^2  +  P2)  —  uo(Po£o  +  po)  =  S(p2^2  —  PqEq)  ■  (59) 


The  state  on  the  left  side  of  the  shock  as  shown  is  the  state  U2.  In  this  case  we  have  po  =  0,  and  this  leads  to  the  conclusion 
that 

u2  =  Uq  =  S,  P2  =  Po-  (60) 

But  if  Ur  is  the  vacuum  state  Uo,  then  p 0  =  0,  hence  P2  =  0.  The  conditions  across  the  contact  discontinuity  show  that 
p  1  =  0.  The  solution  scenario  shown  in  Fig.  2.  If  the  regularity  conditions  (19)  are  imposed,  then  the  zero  pressure  state  is  a 
zero  density  state.  Hence  both  states  Uj  and  U2  are  vacuum  states,  now  connected  to  a  left-going  rarefaction  fan.  Therefore 
the  solution  of  VRP  with  vacuum  right  state  consists  of  left  rarefaction  wave  and  the  tail  of  the  rarefaction  merges  at  the 
vacuum  interface  as  shown  in  Fig.  2(b). 
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Fig.  3.  (a)  Tracking  of  vacuum  interface  boundary  in  the  case  that  the  MVI  moves  across  the  cell  interface  x,+i.  (b)  The  case  of  non-zero  flux  across  the 
grid  interface  x,+1 ,  denoted  by  F1+i . 


Thus  Riemann  problem  for  vacuum  right  state  has  the  left  non-vacuum  constant  states  U i  =  ( Pl,ul,Pl )  and  the  right 
vacuum  state 

U0  =  (p0.K0,po)  =  (0,  0,  0),  (61) 

where  uv  is  the  velocity  of  the  material  vacuum  interface  (MVI).  One  elementary  wave  (isentropic,  rarefaction  fan)  connects 
the  left  non-vacuum  and  the  right  vacuum  state. 


3.2.1  Numerical  flux  approximation 

Numerical  approximation  of  gas  flow  near  the  gas-vacuum  boundary  can  give  rise  to  severe  difficulties  in  the  numerical 
scheme,  associated  with  zero  or  near-zero  density  used  in  mathematical  model  of  the  equation  of  state.  The  Munz  approxi¬ 
mation  scheme  [1]  computes  the  flux  at  the  gas-vacuum  interface  and  simultaneously  tracks  the  vacuum  material  interface. 
The  original  algorithm  was  developed  for  the  case  of  inert,  perfect  gas.  We  use  the  Munz  flux  approximations  and  extend 
the  algorithm  updates  to  multi-component  reacting  flow  uses  non-ideal  EOS  forms  that  are  commonly  used  in  explosive 
modeling.  First,  we  briefly  describe  the  Munz  flux  formula,  for  completeness  and  to  clarify  the  original  presentation. 

Let  x "  be  the  location  of  the  material  vacuum  interface,  with  velocity  u"  at  the  current  time  t  =  tn.  The  location  of  the 
MVI  at  time  tn+ 1  is  estimated  by  the  advance 

x"+1  =xv  +  (Wi  -  tn)u "  .  (62) 


The  numerical  flux  F;+i  at  the  (i  +  T)th  cell  boundary  is  computed  based  on  the  estimated  position  of  the  MVI,  see  Fig.  3(a). 
If  the  MVI  does  not  cross  the  cell  interface  at  Xj+i  at  time  tn+i,  such  that  x"+1  <  Xj+j,  then  the  interface  lies  in  the  vacuum 
during  the  time  step  and  the  flux  F,+1  is  zero.  If  the  MVI  does  cross  the  cell  interface  at  x,+1  at  time  t„+i  (as  shown),  such 
that  x"+1  >  Xj+i  then  the  flux  is  computed  by  means  of  solving  a  local  Vacuum  Riemann  Problem  (VRP)  centered  at  x  =  x". 
The  initial  data  for  the  VRP  (for  a  right  going  MVI)  is 


U(x,  0) 


UL  for  X<0 
U0  for  X>0 


(63) 


where  Uo  is  the  vacuum  state  and  U(.  is  a  non-vacuum  left  state. 

The  VRP  has  the  vacuum  state  on  Uj;  =  Uo  on  the  right.  For  the  left  state  U^,  Munz  proposed  an  average  that  is  weighted 
between  the  cell-centered  states  on  the  left  and  right  of  the  x,  interface 

_  ^fl  _  • 

Ui=aU?+1/2  +  (l  -QOU"  2  with  a  = — ? - L.  (64) 

X;+l  —  xi 


The  cell  between  x,  and  X;+i,  has  vacuum  in  the  interval  [x",Xj+i],  with  value  U"+1/2  in  that  region.  The  average  of  the 
non-vacuum  state  over  the  entire  cell  is  U"+1/2 's  defined  by  U"+1^2(x"  —  x,-)  =  U"+1/2(x;+i  —  x,)  such  that 


Ul  =  U?+1/2  + 


*i+i  -  K 
Xi+l  -  xi 


u 


n 

i— 1/2 


(65) 


is  used  for  the  left  state. 
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Munz’  approximate  solution  to  this  VRP  is 


Ul 

for 

X/r  <  uL  -cL 

U(x,  t)  =  ■ 

U Lfan 

for 

ul  -  cL  <  X/r  <  uv 

(66) 

U0 

for 

X/r  >  uv 

where  all  the  states  are  constants.  The  constant  state  Uyan  is  obtained  by  using  the  integral  form  of  the  conservation  law, 
by  integrating  the  conservation  law  over  the  region  bound  by  X  =  —Ax/2,  X  —  uvr  and  At.  See  Fig.  3(b),  gives 


uvAt 

/ 

-Ax/2 


,  Ax 
U(x,  A t)dX  =  —  UL  +  F(UL)  At . 


(67) 


The  state  U(x,  At)  has  the  value  Uyan  for  x  e  [(14  —  cp)At,  uv  At]. 

The  state  in  that  fan  region  is  approximated  with  an  averaged  value  Uyan.  If  one  applies  the  integral  conservation  law  to 
the  region,  x  e  [—Ax/2,  Xj+i]  and  t  e  [t„,  t„+i]  one  obtaines 


(uL-cL)At 


f 

-Ax/2 


U(x,  A t)dX  + 


uv  At 

/ 


U(x,  A t)dX  =  ^UL  +  F(UL)  At , 


(68) 


(uL-cL)At 


Ax 


Ax 


(uL  -  cL)AtUL  +  — UL  +  [uv  -  (uL  -  cL)]  AtUyan  =  —  UL  +  F(Ul)  At , 

which  obtains  the  formula  for  Uyan 
F(Ui)  -  (ut  -  cL)UL 


U 


Lfan  : 


(69) 


(70) 


Uv  -UL  +  CL 

This  approximation  is  subject  to  the  standard  restriction  that  the  waves  generating  from  an  interface  do  not  interact  within 
the  cell  by  waves  from  adjacent  cell  interfaces,  which  implies  the  CFL  condition 

Ax/2 


At  < 


(71) 


where  Smax  is  the  maximum  wave  velocity,  i.e.  max[  \ul  —  Cl|,  |uv|  ]. 

The  next  step  is  to  use  the  approximate  solution  (66),  (70)  to  calculate  the  flux  at  x, +1,  when  x"+1  >  x;+i,  see  Fig.  3(a). 
The  left  non-vacuum  states  lie  in  the  interval  [—Ax/2,0]  and  the  vacuum  takes  place  in  the  shifted  interval  [0,  Ax/2]. 
Integrating  the  integral  conservation  law  across  the  box  from  [— Ax/2,x!+i  —  x"],  and  [0,  At],  leads  to 


■1  +  1 

/ 


U(x,  A t)dX  =^-UL  +  F(UjA  At  -  Fi+i  At . 


(72) 


—Ax/2 

This  can  be  recast  as 

(u[-c[)At 

U(x,  At)dX  + 

-Ax/2  (u[-c[)At 


L  L 

/ 


X/+ l-*v 

/ 


U(x,  A t)dX  =  +  F(U£)  At  -  Fi+1  At , 


(73) 


>■  _  1-11+1 


U(x,  At)  = 


Uj.  for  —  Ax/2  <  X  <  (ul  —  Cl) At 

Uyan  for  (uL  -cL)At  <  X  <  (xi+1  -x") 

From  the  approximate  solution  (66),  and  Eq.  (72),  we  obtain 
Ax 


Ax 


(uL  -  cL)  AtUj.  +  — Ul  +  [xi+i  -  xnv  -  (uL  -  cL)  AtJUL/^  =  — U L  +  F(UL)  At  -  Fi+i  At . 


(74) 


(75) 


By  separate  consideration  of  the  cases  when  the  left  rarefaction  boundary  is  such  that  u[  —  crL  >  0  or  u[  —  c[  <  0,  the  flux 
F;+i  at  x!+i  computed  from  (75)  is  written  concisely  as 

Fi+1  =  F(U£)  -  -Tmin  [(u[  -  c[)At,  xw  -  x"] 

-  T;max[0,  xi+1  -  xnv  -  ( u[  -  c[)At]  U[fan  . 


(76) 
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In  summary,  the  method  updates  the  interface  location  and  determines  approximate  fan  states  in  the  MVI  adjacent 
cell,  consistent  with  the  integral  conservation  law,  then  updates  the  flux  at  the  vacuum  interface  boundary.  The  tracking 
method  is  only  activated  when  the  density  of  left  (or  right)  local  Riemann  problem  has  a  neighboring  vacuum  state.  If  any 
neighboring  state  is  not  a  vacuum  state,  the  algorithm  computes  the  flux  using  the  general  Riemann  solver  (in  our  case  the 
Roe-solver  for  general  equation  of  state).  A  simple  criterion  is  required  to  establish  a  vacuum  (or  extremely  low  pressure 
state).  In  our  case  a  density  cut  off  is  used,  and  for  the  applications  discussed  in  the  following  section  sets  the  vacuum 
criteria  to  10-9  g/cm3.  When  the  ideal  equation  of  state  is  used  and  no  reaction  occurs,  the  velocity  xv  of  vacuum  interface 
can  be  computed  analytically  as  shown  by  Munz.  However  when  non-ideal  EOS  for  reactive  flow  is  used,  the  computation 
of  vacuum  interface  is  not  simple  and  depends  on  how  the  reaction  quenches  during  the  isentropic  expansion.  Furthermore 
an  expression  for  the  pressure  in  the  isentropic  fan  cannot  be  obtained  in  closed  form  and  numerical  integration  is  required. 
We  describe  these  essential  considerations  in  the  next  section. 


3.3.  Considerations  and  the  limiting  conditions  near  the  vacuum  edge  for  multi-component  reacting  flow 


In  previous  section,  we  described  the  basic  tracking  algorithm.  Here  we  add  the  considerations  required  to  employ  the 
algorithm  for  a  multi-component  reacting  flow  with  mechanical  and  thermal  mixture  EOS  form  expressed  as  of  the  form 
p  =  p(v,  e,  >.;)  and  p  =  p(v,  T, \j).  Menikoff  [22]  worked  the  details  for  a  complete  MG  EOS,  which  we  choose  to  model  the 
component  EOS  forms  for  the  discussion  in  this  section.  We  assume  that 

e,(p,  vO  =  ef(Vj)  +  £-(p-  psi (vt))  +  e0i ,  (77) 

1  i 

where  we  assume  that  pf(v,)  and  e?(v,)  are  specified  reference  functions,  and  in  the  case  when  an  isentrope  is  used  as  a 
reference  curve,  e*  is  such  that  3e?/3vj  =  —  pf (v,).  Both  Tj  and  eo;  are  assumed  constants.  If  the  specific  heat  is  constant 
and  the  reference  curve  is  an  isentrope,  Menikoff  shows  that  the  thermal  equation  of  state  is  given  by 


p  =  pf(Vi)  + 


Cyi  r, 
Vi 


T-T  o 


(78) 


The  volume  vq  is  a  reference  volume.  One  can  also  work  out  the  form  for  the  component  entropies,  to  obtain 


Slip,  Vi) 


s0i  +  Cvi  In 


vo 


k (Cyi  T 0 


[p  -  pf(Vf)] 


(79) 


The  mixture  energy  and  volume  and  entropy  are  given  by 

e  =  ^e;7.;,  v  =  ^v;7.,-  and  s  =  ^s,A.,-. 
i  i  i 


(80) 


In  the  region  near  the  vacuum  edge  the  flow  is  highly  rarefied,  and  nearly  collisionless  so  it  is  reasonable  to  assume  that 
the  reaction  rates  are  zero  (as  is  the  basic  premise  of  the  TOFMS  experiments),  and  hence  isentropic.  The  task  of  working 
out  the  structure  of  the  flow  in  the  fan  is  the  task  of  finding  the  isentrope  in  the  vacuum  fan  solution.  We  briefly  discuss  the 
nature  of  that  calculation  using  these  non-ideal  EOS  forms,  without  further  assumption,  except  that  of  isentropic  expansion. 

The  analysis  of  the  left  rarefaction  fan  at  the  right  edge  of  the  domain  then  requires  an  expression  for  the  p,  v  isentropes 
for  the  mixture.  In  turn  this  requires  finding  expressions  (or  numerically  calculating)  the  values  for  the  partial  volumes 
v,(p,  v)  of  components  in  the  mixture.  The  mixture  entropy  is  then  given  by 


s  =  J2s(Vi(P.v).p)Xi 
i 


(81) 


and  is  found  as  a  function  of  p,  v  and  X, ,  in  general.  The  determination  of  v;(p,  v)  depends  on  the  forms  of  the  mixture 
closure.  If  the  mixture  closure  uses  the  PT-EQB,  an  iterative  numerical  solution  is  required  even  for  a  simplified  version  of 
the  MG  form.  But  that  can  always  be  done  in  principle,  as  the  formulation  is  complete.  The  sound  speed  squared  on  the 
isentrope  is  given  by  c2  =  3p/3p|s,  and  can  be  calculated 


2  2 
c  =  v 


de  , 


P  +  )  / ^Hv,a 


de  . 


3  P 


(82) 


to  obtain  c2(p,  v,  Xf).  The  constant  entropy  in  the  left  state  sp  (say)  can  be  used  in  (81)  to  generate  the  p,  v  isentrope  that 
applies  in  the  rarefaction  of  the  form  sl  =s(p.v,Xi).  In  turn  this  can  be  used  to  construct  the  solution  and  compute  the 
MVI  velocity. 

An  equivalent  and  computationally  simple  procedure  is  available  by  a  directly  integrating  states  through  the  fan  as 
follows.  From  the  consideration  of  the  Generalized  Riemann  Invariants  for  the  system  (46)  for  the  rarefaction  wave  pt\  = 
u  —  c  we  have 
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dp  = 


d(pu)  d(pe) 


(83) 


u  —  c  H  —  uc 

The  first  equality  in  Eq.  (83)  gives  a  formula  for  the  MVI  velocity,  uv  where  Pr  =  0  and  uR  =  uv  as  shown  in  Fig.  2(b). 

PL 


uv  =  uL+  /  -dp. 

J  P 


(84) 


The  second  equality  from  Eq.  (83),  obtains  the  isentropic  change  of  the  internal  energy  in  the  isentropic  expansion  and  is 
written  as 


de  =  —=dp . 


(85) 


This  equation  can  be  used  to  derive  a  differential  equation  for  p(p,Xi).  As  an  example,  consider  the  equation  of  state 
defined  by  using  the  JWL  EOS  forms  for  two  components  as  defined  by  Eq.  (29).  Using  the  form  (34)  the  energy  can  be 
written  as 

e(p,  v ,  X)  =  A  +  G  p  (86) 

where 


vr 


A(vr,  vp,  A.)  =  (1  —  A.)  I  eR(vR) - Pr(vr)  +M  esP(vp) - pP(yp) 


Wr 


vP 


Wp 


and 


Vp  Vr 

G(vr,  vp,  X)  =  X  —  +  (1  —  X)——  . 

COp  Mr 


(87) 


(88) 


The  closure  condition  is  used  to  determine  the  component  specific  volumes  vP  and  vr.  The  PT-EQB  condition  for  the 
thermal  EOS  form  is 

T  =  Tr(vr,  p)  =  Tp(vP,  p) ,  (89) 

where  Tr  and  Tp  are  obtained  from  the  thermal  EOS  for  the  components.  With  <t>  =  vP/vR,  and  Eq.  (30)  leads  to 
v  v<J 


Vr  = 


Vp  = 


(90) 


1-A  +  A4>’  '  1-A  +  A.4> 

Eq.  (89)  determines  a  value  of  4>(p,  v)  for  fixed  X.  Given  a  value  X,  the  functions  A  and  G  are  functions  of  v  and  p. 
Then  if  we  substitute  Eq.  (86)  in  (85)  one  obtains  a  differential  equation  for  the  pressure  as  a  function  of  density  for  fixed 
composition,  that  holds  in  the  rarefaction  wave 

dp  1 

^  =  c+34+P£ 


'dp 


rr  j__ /ac\l 

dA~ 

.Ip2  \9  p) 

p~Jp_ 

subject  to  P(Pl)  =  Pi¬ 


rn 


Note  that  the  c2  =  dp /dp  on  the  isentrope.  The  interface  velocity  is  then  computed  from  (84)  as 


U!  =  UL+  /  - 
P 


PL 

f~ 

J  P 


r  I  3/1  I  n  3C 

G  +  3p  +P3F 


3  P 


1  /  9G 

P2  \3P 


3A 

If, 


dp. 


(92) 


The  pressure  p  in  the  integrand  of  Eq.  (92)  is  the  solution  of  the  ODE  (91). 

Generally  the  ODE  for  p  and  the  integral  for  uv  must  be  computed  numerically.  However  if  the  simpler  P-EQB  closure 
is  used  to  compute  the  component  volumes  as  in  Stewart  et  al.  in  23],  then  <I>  does  not  depend  on  pressure  and  the  ODE 
(91)  simplifies  to 


dp  1 

1 

/  dG\ 

dA~ 

dp  G 

p2 

\dp)\ 

[P_^_ 

subject  to  p(Pl)  =  Pl, 


and  can  be  integrated  analytically  to  obtain  the  solution  for  p 

P(P)  =  Pi  +  r pr+1  (a(Px)p~t  -  A{p)p-r  +  r  J  Ap-1'-1 


dp 


(93) 


(94) 


where 


T  =  (l  -A  +  A4>)/ 


Mr, 


Mr 


so  that  G  =  p. 


(95) 


S.  Choi  et  al. /Journal  of  Computational  Physics  296  (2015)  158-183 


171 


At  the  limit  of  zero  density,  one  must  check  the  Liu-Smoller  conditions  (19)  to  make  sure  that  the  pressure  and  sound 
speed  vanish  as  the  density  vanishes.  The  solution  satisfies  the  end  condition  p(0)  =  0.  The  derivative  is  given  by 


dp 

/  =  (r  +  D 

dp 


PL 


Pl 


r+i 


+  r^<Pt) +r2 
Pi 


Pl 

fi 


>r+i 


,  \  r  dA 

dp  \Pr-Tp—  -nr  +  iM(p), 


(96) 


and  the  solution  satisfies  the  Liu-Smoller  conditions  (19)  if  A  =  o(p 1  +1)  and  A!  —  o(p'  +1)  as  p  — >  0. 

Physical  consideration  of  the  near  collision  limit  of  a  vacuum  would  admit  the  reasonable  approximation  that  each  EOS 
for  each  gaseous  component  of  the  mixture  limits  to  an  ideal  gas  form.  For  each  component  in  the  limit  of  infinite  volume, 
or  zero  density  we  would  require  that 


cvj  r  j 


Vi 


To(  — 
V  Vi 


0  as  pj  =  1  /  Vf 


0. 


so  that  in  the  limit  of  low  densities  the  ideal  EOS  forms  apply 
PV,-  ,  .  RJ 


e,(p ,  v,) • 


+  e0i ,  and  p  ~ - , 


(97) 


(98) 


(Yi-  1) 

where  — »•  R,-  (the  ideal  gas  constant  weighted  by  the  molecular  weight  of  the  component)  and  T,  — >  Yi  ~  1-  One  can 

view  these  additional  requirements  as  a  modeling  statement  that  reflects  the  correct  physics  or  simply  a  physically-based 
regularization  for  the  model  that  is  used  only  near  the  MVI.  In  either  case  the  additional  requirement  that  each  gaseous 
component  EOS  form  limits  to  ideal  EOS  form  when  PT-EQB  is  used  as  the  mixture  closure  model,  will  guarantee  that  the 
Liu-Smoller  conditions  are  met  which  we  show  next. 

For  the  conditions  near  the  MVI,  the  state  will  be  such  that  the  EOS  forms  of  the  components  of  the  mixture  limit  to 
the  ideal  EOS  forms.  For  fixed  pressure  and  temperature,  we  can  take  the  expression  for  the  partial  volumes  and  sum  them 
over  the  components  to  obtain 


T  RT 
P 

where  R  —  is  the  mass  weight  ideal  gas  constant.  Then  specific  energy  of  the  multi-component  mixture  limits  to 


(99) 


PV 


Y 


Y  +  e‘°^'  ’ 


with 


1 


T-? 


Ai  Ri 


(100) 


and  where  we  have  used  T  ■ 
sound  speed  squared 


Y  ~  1  (Yi  ~  1)  R 

pv/R  and  identified,  and  y  is  a  function  of  the  composition.  Then  the  standard  formula  for 


c  =  y  pv 

and  the  isentrope  in  the  fan  expansion  is  given  by 


(101) 


(102) 


where  the  constant  k  is  also  a  function  of  the  frozen  composition,  which  determines  y.  If  these  limiting  forms  for  the 
EOS  are  used  then  the  Liu-Smoller  conditions  are  satisfied  since  p/py  —  k  —  Pl/(p[)  >  0  and  c2  =  yp/p  implies  c  = 
<Jy  k  p(v-1)/2  with  the  properties  that 

c{p)  >  0  for  p  >  0  and  c(0)  =  0,  (103) 


and 


Pi 

/: 


cdp=2^Pri)/2 

p  y- 1 


bounded . 


(104) 


It  is  important  to  compute  the  integral  (92)  carefully,  since  there  is  often  an  abrupt  change  in  the  states  through  the 
rarefaction,  especially  when  a  shock  wave  first  hits  the  MVI.  To  avoid  singularity  occurring  from  lower  limit  of  the  integrand 
in  (92),  we  always  split  the  integration  into  two  pieces  as  shown  in  Fig.  4.  We  define  a  cut-off  density  e  that  is  small 
enough  to  ensure  that  the  ideal  gas  limiting  form  is  a  good  physical  approximation  to  the  flow,  and  write  the  integral 


Pl  €  Pl 

(  -dp=  f-dp+  f1 
J  P  J  P  J  i 


-dp . 


(105) 


172 


S.  Choi  et  al.  /  Journal  of  Computational  Physics  296  (2015)  158-183 


Fig.  4.  Two  ranges  of  integration  of  Eq.  (105). 


For  a  given  specified  value  of  e,  with  0  <  e  <<  1,  the  first  term  of  (105)  is  evaluated  with  Eq.  (104) 


[  c_dp  =  2Jlltf- m 
J  p  y- 1 

0 


(106) 


and  the  second  integral  is  computed  numerically  using  the  Romberg’s  method  [24], 

As  an  example,  the  JWL  EOS  forms  are  MG  forms  that  are  commonly  used  to  model  explosives  and  they  have  the 
property  that  they  limit  to  ideal  EOS  forms  as  the  component  densities  vanish.  A  JWL  EOS  for  a  component  takes  the  form 

ps.  =  Aie~R i‘(v'/vo)  +  gie-R2i(vi/v0)  +  RiTo  { X°Y ' 

'  '  V0  V  Vi  ) 


and  as  p,-  =  1/v,  ->  0,  it  follows  that  (97)  holds.  Thus  the  JWL  EOS  behaves  like  the  ideal  gas  near  vacuum  region  or  density 
is  very  small,  which  is  a  feature  built  into  the  JWL  EOS  calibration. 

If  the  non-ideal  EOS  forms  for  the  components  do  limit  to  the  ideal  EOS  forms  in  the  limit  of  vanishingly  small  density, 
then  the  classic  fan  analysis  will  hold  at  the  extreme  edge  of  the  fan.  Even  then  we  make  the  observation  that  the  computed 
MVI  velocity  can  be  very  sensitive  to  the  limiting  values  of  the  flow  near  the  fan  when  the  reaction  ceases,  since  the 
composition  of  the  frozen/quenched  products  determines  the  mass  fraction  weighted  value  of  y  of  the  mixture. 

Close  to  the  MVI  assume  that  the  expansion  to  vacuum  occurs  by  a  left-going  isentropic  fan  seeded  by  a  constant  left 
state  that  lies  in  a  region  a  nonzero  density,  but  where  the  ideal  EOS  forms  are  applicable.  This  would  be  a  point  close  to 
the  MVI  edge,  but  to  the  left  of  it  (say).  Then  the  particle  velocity  and  speed  relation  hold  between  the  defined  left  state 
and  any  state  to  the  right  toward  the  MVI,  i.e. 

2c  2  Cl 

u  + - -  =  uL  + - Y-  (108) 

y  —  1  y  —  1 

Eq.  (108)  can  be  used  to  illustrate  why  the  tracking  of  vacuum  interface  numerically  can  be  very  difficult  if  done  without 
knowledge  of  the  analytical  character  of  the  fan  solution. 

Since  the  fan  boundary  at  the  MVI  is  characteristic,  then  uv  =  X/r,  and  this  is  consistent  with  setting  u  —  uv  and 
p  =  c  =  0  at  the  vacuum  interface.  Then  (108)  becomes 

uv  =  X/r  =  uL  +  (109) 

Eq.  (109)  could  be  used  to  compute  uv  near  the  interface,  but  since  the  Liu-Smoller  requirement  insists  that  in  a  (compu¬ 
tationally)  small  neighbor  of  the  interface  cl  ->  0,  under  resolution,  then  this  relation  is  simply  the  identity  uv  —  X/ r  =  ui. 

But  now  consider  Eq.  (108)  in  the  same  small  neighborhood,  and  consider  another  neighboring  state  to  the  right  of  the 
same  left  state  in  the  left  fan.  We  identify  this  state  with  a  vacuum  cut-off  density  d  (say)  to  designate  being  close  enough 
to  the  MVI.  Replace  uv  =  X/r,  as  before,  but  estimate  the  sound  speed  which  is  a  function  of  p  at  the  vacuum  cut  off  as 

/  d  Vy~1)/2 

c(d)  =  cL^—J  (110) 

in  which  case  (108)  becomes 


es  _  Yie-RiiOVv0)  ,  lLe-R 2i(v,-/v0)  ,  R‘T°  ( ^ 

'  Ru  Rii  (Yi  -  1)  V  v/ 


(n- 1) 


(107) 


uv  =  uL  + 


2  CL 


2  CL 


y  -  l  y  -  l 


^yk-D/2 

Pi) 


(in) 


There  are  clearly  many  ways  to  go  wrong  if  the  local  approximations  near  the  vacuum  interface  are  not  accurate,  incon¬ 
sistent  with  the  Liu-Smoller  conditions  or  inconsistent  with  the  fan  structure.  The  greatest  sensitivity  in  the  formulas  are 
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Fig. 5.  The  plot  of  interface  velocity  function  uv(d)  shown  in  Eq.  (Ill)  with  pL  =  3.0e— 8  [g/cm3],  pL  =  3.6e— 8  [GPa],  uL  =  16.5  [mm/ps], 

illustrated  by  the  dependence  on  the  local  value  of  y.  As  y  — »•  1  the  interface  velocity  function  (111)  becomes  singular,  and 
thus  is  very  sensitive  to  y.  Fig.  5  shows  the  plot  of  velocity  uv(d )  of  Eq.  (Ill)  as  a  function  of  density  d  for  different  values 
of  y;  i.e.  y  =  1.2, 1.4, 1.8,  and  2.0.  The  left  state  variable  pi,  pi  and  Ui  shown  in  Eq.  (Ill)  was  set  to  actual  computational 
values  at  time  t=  1.5  ps  used  in  the  PETN  problem  in  Section  4.2.  As  shown  in  the  figure,  the  interface  velocity  increases 
rapidly  and  the  slope  of  uv(d)  becomes  very  large  as  y  decreases  toward  one.  Also  it  is  clear  from  the  dependence  of 
pi  in  the  denominator  of  Eq.  (Ill)  that  if  the  monotonically  decreasing  ratio  of  (d/ pi)  is  not  maintained  correctly  in  the 
approximation  scheme,  even  for  fixed  cp  that  the  prediction  of  the  interface  velocity  can  be  made  highly  inaccurate.  Since  y 
is  a  function  of  the  mass  fraction  as  seen  in  Eq.  (100),  it  is  dependent  on  the  composition  of  the  products  near  the  vacuum 
edge.  Thus  the  details  of  the  kinetic  rate  law  and  how  the  reactions  quench,  does  matter  in  the  calculation  of  the  MVI 
velocity  uv. 

The  conclusions  from  the  simple  considerations  above  are  independent  of  the  grid  resolution.  The  vacuum  tracking 
algorithm  uses  the  fan  solution  that  maintains  the  correct  limiting  forms,  derived  from  the  EOS  forms  of  the  underlying 
model  as  the  vacuum  state  is  approached.  The  requirement  that  the  EOS  forms  limit  to  the  ideal  EOS  forms  and  that  the 
reaction  stops,  is  easy  to  implement  and  is  physically  consistent.  The  JWL  EOS  form  is  an  example  of  a  non-ideal  form  that 
has  the  correct  properties  and  the  Liu-Smoller  condition  are  automatically  satisfied.  Any  vacuum  tracking  scheme  that  would 
work  similarly,  must  employ  consideration  of  the  EOS  forms  used  for  the  mixture  in  the  vicinity  of  the  vacuum,  otherwise 
large  and  likely  catastrophic  inaccuracies  will  occur  that  will  destroy  any  simulation  that  has  a  vacuum  or  extremely  low 
pressure  regions.  In  the  next  section,  we  demonstrate  application  of  the  algorithms  with  examples  that  explicitly  uses  ideal 
and  nonideal  EOS  forms  for  the  mixture  of  reactants  and  their  products. 

4.  Examples  and  applications 

In  this  section,  we  demonstrate  the  applicability  of  our  algorithm  through  two  representative  numerical  simulations.  The 
first  example  uses  the  ideal  equation  of  state.  The  second  example  is  representative  of  the  TOFMS  experiments  which  uses 
the  non-ideal  JWL  EOS  forms  for  reactants  and  products. 

4.1.  Example  1:  An  explosive  modeled  by  the  ideal  gas  EOS  form  whose  products  expand  into  a  vacuum 

In  many  theoretical  studies  of  high  explosives  (that  are  condensed  solids  at  their  initial  state),  the  flows  are  modeled 
with  the  ideal  EOS  forms  for  a  mixture  of  reactants  and  products.  The  reader  is  referred  to  Fickett  and  Davis  [12].  We  used 
values  listed  in  their  text  for  this  example  [12,  Chapter  2,  C2,  p.  46],  The  EOS  is  given  by 

e(p,v,X)  =  -^--QX.  (112) 

7-1 

The  density  in  the  unreacted  explosive  is  assumed  to  be  po  =  1.6  g/cm3.  For  the  choices  y  =3.0,  and  Q  =4.5156  MJ/kg, 
the  Chapman-Jouguet  (CJ)  detonation  velocity  and  pressure  are  Dq  =  8.5  mm/ps  and  pcj  =  28.9  GPa,  respectively. 

The  computational  domain  is  taken  to  be  50  mm  long.  The  boundary  at  x  =  0  is  an  outflow  boundary  that  enforces 
zero  gradients.  The  material  in  the  region  0  <  x  <  20  mm  is  assumed  to  be  a  uniform,  completely  reacted,  CJ  detonation 
state  with  the  state  vector  as  Uq  =  (p.u.p.A)  =  (pq,  uq,  pq,  1).  In  the  strong  shock  approximation,  the  CJ  states  are 
pq  =  (y  +  \)/y)po,  uc  —  Dq/(y  +  1),  and  pq  =  poD^/(y  +  1).  The  material  between  20  <  x  <  30  mm  is  assumed  to  be 
at  a  uniform,  completely  unreacted,  motionless  state  at  an  initial  density  po,  but  at  zero  pressure  with  Uhe  =  (p,  u,  p,  A.)  = 
(po,  0,  0,  0).  The  region  to  the  right  of  x  =  30  mm  is  assumed  to  be  the  vacuum  state  Uo  =  (0.  0,  0,  0).  Thus  the  initial  values 
are  given  by 
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Table  1 

Initial  conditions  for  the  simulation  in  Example  1. 


State 

Reacted  CJ  state 

Unreacted  CJ  state 

Vacuum  state 

U 

0  <  x  <  20 

20  <  x  <  30 

20  <  x  <  30 

p 

2.133  g/cm3 

1.6  g/cm3 

0 

u 

2.125  mm/ps 

0 

0 

p 

28.9  GPa 

0 

0 

U(x)  = 


Uq  (reacted  CJ)  if  x  <  20 
Uhe  (unreacted  HE)  if  20  <  x  <  30 
Uq  (vacuum)  x  >  30 


and  listed  in  Table  1.  The  reaction  rate  in  the  explosive  is  taken  as 


(113) 


A  =  (5.0)('l -A)05  (— ^  ps~\  (114) 

VPcj/ 

The  initial  high  pressure,  density,  particle  velocity  state,  Uq  introduces  a  shock  into  the  unreacted  HE  at  the  x  =  20  mm 
interface  and  initiates  a  detonation  wave,  that  builds  up  and  runs  in  the  unreacted  HE  towards  the  vacuum  interface  initially 
at  x  =  30  mm.  Since  the  pressure  in  the  unreacted  HE  is  zero,  there  is  no  initial  expansion  into  the  vacuum,  so  the  MVI  at 
x  =  30  mm  is  motionless  until  the  detonation  shock  hits  the  interface  and  the  detonation  products  start  to  vent  into  the 
vacuum  region.  The  initiated  detonation  becomes  nearly  steady,  with  a  von  Neumann  spike  shock  pressure  of  about  55  GPa 
pressure  at  approximately  t  =  1.3  ps.  At  that  time  it  encounters  the  previously  static  vacuum  interface.  After  the  detonation 
shock  and  the  following  flow  in  the  detonation  reaction  zone  hits  the  vacuum  interface  at  x  —  30  mm,  the  detonation 
shock  disappears  and  a  rarefaction  wave  is  sent  backwards  to  the  left  into  the  gas  flow.  Since  reaction  rate  in  the  explosive 
is  a  function  of  the  pressure,  the  reaction  rate  in  the  fan  region  drops  as  the  pressure  drops.  Figs.  6  and  7  show  the  state 
variables  at  different  times  and  show  the  formation  of  a  detonation,  followed  by  the  expansion  of  products  as  they  vent  into 
the  vacuum  region.  The  simulation  uses  5000  grid  points  for  the  entire  domain,  and  exhibits  cleanly  defined,  well-resolved 
features  and  a  sharply  defined  MVI  motion. 


4.2.  Example  2:  detonation  of  pentaerythritol  tetranitrate  (PETN)  whose  products  vent  to  vacuum 


Next  we  describe  a  simulation  of  the  expansion  of  PETN  products  into  a  vacuum,  that  represents  a  facsimile  of  TOFMS 
experiments.  The  JWL  EOS  parameters  for  reactants  and  products  are  chosen  to  match  experimental  data  obtained  from 
macroscopic  detonation  experiments  on  PETN  [26],  These  experiments  include  those  that  determine:  a)  the  unreacted  shock 
Hugoniot  or  Up-Us  relation,  b)  the  overdriven  detonation  shock  Hugoniot  for  the  products,  and  c)  the  products  expansion 
experiments  that  determine  constants  for  the  products  JWL  EOS  [27,25,28,26  .  The  reaction  rate  parameters  are  determined 
from  shock  initiation  experiments  that  determine  the  distance  and  time  to  detonation  as  a  function  of  the  shock  input 

pressure  into  the  sample.  The  reader  can  find  details  on  how  the  equation  of  state  and  rate  law  parameters  are  fit  to 

explosives  in  recent  review  papers  [29, 30], 

We  use  model  parameters  reported  by  Tarver  et  al.  in  [26],  The  calibrated  parameters  of  the  JWL  EOS  forms  for  the 
reactant  and  product  are  given  in  Table  2.  The  reference  density  is  given  by  p o  =  1.778  g/cm3.  The  JWL  product  equation 
of  state  parameters  determines  the  CJ  states  which  are  found  to  be 

Dq  =  8.23  mm/ps,  vq  =  0.42  [cm3/g],  pq  =  31.5  GPa  (115) 

The  distance  to  detonation  for  full  density  PETN  is  reported  as 

l°gio  p  =  ci  —  b  log^  x*  (116) 

where  p  [GPa]  is  the  shock  input  pressure  and  x*  [mm]  is  observed  distance  to  detonation  from  the  explosive  edge,  and 
parameters  a=  1.1767  [GPa],  b  =  0.2432.  We  used  a  simple,  pressure  dependent  reaction  rate  model  with  rate  parameters 
chosen  so  that  our  PETN  model  reproduces  the  reported  the  run  to  detonation  data,  (116).  The  rate  is  given  by 

^  =  80  (1  —  7.)0'5  ( ~^\  ps^1.  (117) 

dt  \pc j) 

Fig.  1  shows  a  schematic  of  the  setup  for  a  simulation  of  motionless  PETN  that  has  vacuum  on  the  right  side,  and  is 
struck  by  a  impact  piston  on  the  right  side.  Initial  states  (p.p.X)  at  initial  time  t  =  0  are  divided  into  three  regions  and 
are  shown  in  Table  3.  Through  impact  on  the  left  side,  a  shock  is  introduced  into  the  unreacted  material.  If  the  shock  is 
strong  enough  and  the  sample  is  thick  enough,  then  a  steady  state  detonation  forms  in  the  PETN  sample  in  the  region 
X;  <  x  <  xv  before  the  detonation  shock  strikes  the  vacuum  interface,  after  which  the  venting  at  the  interface  occurs.  In 
PETN,  the  steady  CJ  detonation  structure  has  a  reaction  complements  of  approximately  99%  in  0.02  mm,  and  the  final  1% 
of  the  energy  that  relaxes  to  the  CJ  state,  is  released  in  a  reaction  zone  tail  that  is  approximately  0.2  mm.  Therefore  the 


S.  Choi  et  al.  /  Journal  of  Computational  Physics  296  (2015)  158-183 


175 


Fig.  6.  Sequence  showing  the  detonation  hitting  the  MVI  and  the  expansion  of  products  into  the  vacuum  region:  (a)  density,  (b)  pressure,  (c)  velocity, 
(d)  reaction  progress. 
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Fig.  7.  Space  time  contour  plots  that  show  the  same  sequence  as  Fig.  6  showing  the  detonation  hitting  the  MVI  and  the  expansion  of  products  into  the 
vacuum  region:  (a)  density,  (b)  pressure,  (c)  velocity,  (d)  reaction  progress  field  in  space  and  time  domain. 


Table  2 

The  calibrated  parameters  of  JWL  EOS  of  PETN. 


A,  [GPa] 

B,  [GPa] 

Ci 

Wj 

Ri, 

R2i 

eoi  [KJ/g] 

Product  JWL  (i  =  R) 

1032.158 

90.57 

0.0 

0.57 

6 

2.6 

6.074 

Reactant  JWL  (i  =  P) 

1280.0 

-27.058 

0.0 

0.6 

6 

2 

0.0 

Table  3 

Initial  data  for  impact  simulation  Example  2. 


Material 

P  [Gpa] 

p  [g/cm3] 

u  [mm/ps] 

Range 

Piston 

0.0 

2.14 

Uin 

X<Xi 

PETN 

0.0 

1.778 

0.0 

Xi  <  x  <  xv 

Vacuum 

0 

0 

0.0 

x>xv 

computation  requires  high  resolution,  and  the  spatial  grid  size  must  be  no  more  than  one  micrometer  with  approximately 
20  points  in  the  rapidly  changing  detonation  shock  structure.  Note  that  a  multi-material  solver,  described  in  [31],  was  used 
to  properly  model  the  flyer  impact  (of  a  different  material,  specifically  a  plastic  KEL-F)  and  transmission  of  the  impact  shock 
into  the  PETN  sample.  But  the  interior  algorithms  and  the  algorithms  near  the  MVI  are  identical  to  what  is  discussed  in  this 
paper. 

Experimental  detection  of  product  species  flying  through  vacuum  depends  on  many  factors,  that  include  the  thickness 
of  sample  and  how  the  sample  is  shock  initiated.  A  premise  of  the  experiment  is  that  a  steady  detonation  is  established 
by  the  time  the  detonation  shock  hits  the  vacuum  interface  and  the  reactants  are  quenched  by  the  subsequent  rarefaction. 
When  the  PETN  sample  is  initiated  with  an  input  shock  strength  of  5  GPa  in  pressure,  matched  to  the  experimental  pop-plot 
shows  a  steady  state  detonation  is  established  in  about  0.16  ps  at  around  0.5  mm  away  from  the  edge  of  the  sample.  In  this 
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Fig.  8.  The  location  of  impactor,  PETN  sample,  vacuum  and  seven  probe  locations  of  measurement  for  simulation  in  Example  2. 

scenario,  any  sample  thicker  than  0.5  mm  is  sufficiently  thick  to  observe  steady  state  detonation.  Once  the  steady  state 
detonation  on  PETN  hits  the  vacuum  interface  and  the  vented  products  expand  at  a  very  high  speed,  faster  than  16  mm/ps. 

Next  we  discuss  a  set  of  cases  in  which  the  thickness  of  the  sample,  i.e.  xv  —  x;  and  the  impactor  velocity  Uimpactor  (i.e. 
related  to  the  input  shock  pressure)  are  varied.  In  each  case  the  density  p,  pressure  p  and  reaction  progress  X  are  sampled 
at  a  designated  point  far  away  from  the  original  vacuum  interface  locations  and  the  time  histories  are  recorded,  similar  to 
what  would  be  recorded  at  a  sampling  probe  at  a  fixed  point  in  space  in  the  TOFMS  experiment.  The  three  different  cases 
are 


•  Case  I:  xv  —  x;  =  3.0  mm,  Uimpactor  =  2.963  mm/ps 

•  Case  II:  xv  —  X;  =  1.0  mm,  u,mpartor  =  2.963  mm/ps 

•  Case  HI:  xv  —  x\  —  3.0  mm,  U[mpactor  =  2.563  mm/ps 

Fig.  8  shows  the  locations  of  the  impactor,  PETN  (HE  sample)  and  vacuum  and  probe  locations  where  measurements  are 
recorded.  The  probe  locations  are  set  to  x  =  100,  200,  300, 400,  500,  600  [mmj  and  x  =  680  [mm]  measured  from  the  initial 
vacuum  interface. 

Figs.  9  and  10  plot  the  state  variables  in  the  region  0  <  x  <  10  mm.  They  show  the  formation  of  a  detonation,  followed 
by  the  initial  expansion  of  products  as  they  start  to  vent  into  the  vacuum  interface.  Fig.  9  shows  the  solution  profiles  of  the 
state  variables  for  Case  I  at  the  indicated  times  from  t  =  0  to  t  =  0.69  ps.  The  solution  profiles  at  t  =  0.39  ps  corresponds 
to  a  nominally  steady  CJ  detonation  wave  that  is  established  in  the  PETN  sample,  just  prior  to  the  detonation  shock  hitting 
the  vacuum  interface  xv  that  is  initially  located  at  4  mm.  Once  the  detonation  shock  wave  hits  the  vacuum  interface,  the 
products  expand  into  the  vacuum  via  the  action  of  a  rarefaction  wave.  Fig.  10  shows  the  corresponding  contours  of  the 
state  variables  plotted  in  the  x-t  (space  and  time)  for  Case  I.  Figs.  9(d)  and  10(d),  clearly  show  the  region  of  interest 
for  the  experimental  study  where  quenching  of  the  reaction  occurs  and  that  would  be  the  region  where  reactions  of  the 
intermediate  species  would  likely  cease,  since  the  pressure  in  the  rarefaction  fan  region  drops  rapidly. 

Next  we  discuss  the  observations  in  the  full  length  of  the  vacuum  region  x  <  680  mm.  After  the  sample  is  detonated, 
the  products  expand  toward  the  vacuum  and  spread  evenly  in  the  entire  vacuum  region  of  the  experiment.  Fig.  11(a)  shows 
the  density  profiles  in  vacuum  expanding  process  from  times  t  =  5  ps  to  t  —  270  ps  for  Case  I.  The  density  was  monitored 
at  various  locations.  Fig.  11(b)  shows  the  monitored  density  at  the  end  of  the  vacuum  region  at  x  =  680  mm  for  Cases  1, 
II  and  III.  The  density  and  time  of  arrival  of  the  density  is  different  for  all  three  cases,  and  is  affected  by  different  sample 
thicknesses  and  impact  velocities  associated  with  each  case.  Fig.  12  shows  the  pressure  and  particle  velocity  plotted  in  the 
x-t  (space,  time)  plane  for  the  entire  vacuum  regions,  for  Case  I.  The  velocity  of  vacuum  interface  can  be  estimated  to  be 
about  19  mm/ps. 

For  all  three  cases  at  each  probe  location,  we  recorded  the  time  of  arrival  (TOA)  that  corresponds  to  a  time  when 
specified  density  value  was  first  sensed.  The  densities  ranged  from  a  low  monitoring  value  of  10-8  g/cm3,  to  a  maximum 
of  10-3  g/cm3.  The  probes  first  senses  the  low  values  and  later  sense  the  higher  values  and  the  results  for  all  three  cases 
are  shown  in  Fig.  13.  Case  II  has  a  smaller  explosive  sample  thickness  than  Case  I,  but  is  subjected  to  the  same  impact 
velocity.  In  Case  II,  the  detonation  wave  is  not  a  steady  state,  CJ  detonation  when  it  hits  the  initial  vacuum  interface. 
Fig.  11(b)  shows  that  a  smaller  peak  density  is  attained  and  that  TOA  at  the  monitoring  locations  has  larger  values.  Thus  the 
recorded  velocity  of  a  constant  density  front,  measured  by  sampling  the  products  expansion  in  the  vacuum  region  is  slower 
for  Case  II  than  for  Case  1.  Case  HI  has  the  same  sample  thickness  of  the  PETN  sample  as  Case  I,  but  the  sample  is  subject  to 
a  lower  impact  velocity.  In  Case  III,  the  sample  is  thick  enough  so  that  a  steady  state  CJ  occurs,  and  the  peak  pressure  and 
density  observed  at  the  probes  is  quite  similar  to  Case  I,  but  very  slight  shift  of  TOA  of  fixed  density  value  at  the  probes  is 
observed. 

These  results  can  be  summarized  by  computing  the  phase  velocity  of  a  constant  density  front  as  the  products  encounter 
each  probe.  The  results  are  listed  in  Table  4.  The  table  shows  that  the  velocity  of  the  front,  for  a  given  monitored  density, 
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Fig.  9.  Time  sequences  of  detonation  and  expansion  processes  for  Case  I:  (a)  density  [g/cm3],  (b)  pressure  [GPa],  (c)  velocity  [mm/ps],  (d)  reaction  progress. 
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Fig.  10.  (a)  Density,  (b)  pressure  (log10(p)),  (c)  velocity,  (d)  reaction  progress  field  in  space  and  time  domain  for  Case  I. 


Fig.  11.  (a)  Density  profile  expanding  toward  the  vacuum  region,  (b)  monitoring  density  at  x  =  680  mm  for  the  Cases  I— III. 


ranges  from  8  to  19  mm/ps  and  is  a  function  sample  thickness  and  impact  velocity.  Experiments  by  Thomas  et  al.  [32] 
showed  that  the  detonation  products  expanding  toward  vacuum  achieve  terminal  velocities  in  the  range  of  8  to  12  mm/ps. 
Fajardo  et  al.  [33],  reported  the  leading-edge  velocities  of  species  through  the  vacuum  approximately  ««  10  mm/ps  in  their 
experiments  on  aluminum  ablation.  Lundborg  [34]  measured  the  velocity  of  the  shock  front  of  TNT  evacuating  from  cham¬ 
bers  and  the  velocity  increased  from  10  mm/ps  to  20  mm/ps  with  decreasing  pressure.  So  our  simulations  easily  fall  in 
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(a)  x  [mm]  (b)  x  [mm] 


Fig.  12.  Pressure  (a)  and  particle  velocity  (b)  fields  for  Case  I  in  x-t  domain. 

the  range  reported  in  related  experiments,  and  is  a  reasonable  check  of  the  applicability  of  these  algorithms  and  simulation 
techniques  for  the  TOFMS  experiments. 

4.3.  Comparison  of  two  different  closure  models  that  are  used  to  define  the  mixture  EOS 

The  vacuum  tracking  interface  algorithms  described  in  this  paper  can  be  generally  used  for  simulating  the  propagation 
of  products  from  the  reaction  of  multi-component  mixtures  and  can  be  modeled  with  any  EOS  that  obeys  the  Liu-Smoller 
conditions  for  the  mixture.  However  when  a  multi-component  formulation  is  used  for  the  mixture,  one  must  provide  some 
closure  conditions  that  are  in  essence  the  statements  of  local  equilibrium  between  components  within  the  mixture.  For  a 
Gibbs  formulation  one  posits  a  single  pressure  and  temperature  (p.  T)  for  the  mixture,  denoted  as  (PT-EQB).  This  means  that 
the  component  EOS  forms  must  also  share  the  same  mixture  pressure  and  temperature.  But  one  can  close  the  mathematical 
formulation  by  proposing  slightly  different  statements.  A  common  one  is  that  while  there  a  common  pressure  (P-EQB),  the 
ratio  of  the  volumes  of  components  is  either  constant  or  is  a  function  of  the  composition,  in  which  case  the  component 
volume  are  not  pressure  dependent. 

This  model  (P-EQB)  enforces  only  pressure  equilibrium  and  specifies  the  volume  component  ratios  as  a  closure  condition. 
It  is  attractive  since  it  uses  EOS  forms  that  are  linear  in  the  pressure  and  the  EOS  forms  are  closed  explicitly.  This  type  of 
closure  eliminates  the  need  for  an  iterative  cycle  to  determine  the  temperature,  and  might  be  considered  if  one  has  two 
(or  multi-)  temperature  model  for  reactants  and  products  that  are  in  pressure  equilibrium  but  not  temperature  equilibrium. 
However  our  purpose  here  is  simply  to  illustrate  that  changing  the  mixture  closure  condition  does  have  an  effect  on  the 
simulated  results,  since  closure  really  represents  a  different  model  of  the  mixture  EOS.  We  note  that  as  we  changed  the 
closure  models,  the  qualitative  features  of  the  simulated  flow  are  quite  similar  but  the  computed  values,  that  are  important 
when  interpreting  chemical  measurements  made  at  the  probes,  do  change. 

To  examine  the  dependencies  of  the  simulation  results  on  the  mixture  closure  models,  we  implemented  the  pressure 
equilibrium  (P-EQB)  closure  for  the  mixture  equation  of  state  in  our  PETN  model  and  carried  out  the  simulations  for  Cases  I 
and  II.  Comparisons  of  simulated  results  for  Cases  I  and  II,  that  use  the  two  closures  (PT-EQB)  and  (P-EQB)  are  shown  in 
Figs.  14  and  15.  The  TOA  in  Fig.  14  is  plotted  versus  the  probe  locations.  The  detection  velocities  of  the  fixed  density  of 
10-s  g/cm3  for  Case  I  are  listed  in  Table  5.  The  velocities  for  the  two  closures,  range  respectively  from  8  to  19  mm/ps  and 
9  to  13.5  mm/ps. 

Fig.  15  shows  a  comparison  of  the  density  recorded  at  x  =  680  mm  for  the  two  closure  models.  Both  models  show  the 
same  tendencies  and  approximately  the  same  densities,  however  there  are  noticeable  differences  in  the  time  shift  and  peak 
density  at  a  fixed  monitoring  point.  The  pressure  equilibrium  (P-EQB)  model  shows  higher  peak  density  but  slower  time 
of  arrival  for  Cases  I  and  II.  These  differences  are  entirely  attributed  to  the  differences  in  the  closure  models.  Thus  careful 
physical  chemistry-based  considerations  must  be  taken  into  account  to  make  accurate  quantitative  predictions  for  this  class 
of  simulations. 

5.  Conclusions 

We  have  illustrated  how  the  original  Munz’s  vacuum  tracking  algorithm  for  a  single  ideal  gas  can  be  extended  to  a  gen¬ 
eral  multi-component  reacting  flow,  with  an  Euler  solver  method  that  uses  a  linearized  Riemann  solver  for  a  Mie-Gruneisen 
EOS  form.  Our  implementation  uses  a  fairly  standard,  cell  based  TVD  scheme  with  source  splitting  to  solve  the  reacting 
Euler  equations.  We  have  shown  that  for  a  non-ideal  EOS  model,  our  method  can  be  implemented  with  different  mixture 
closure  models  for  pressure  and  temperature  equilibrium,  or  just  pressure  equilibrium  that  uses  an  auxiliary  condition  such 
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(b)  monitoring  position  [mm] 


Fig.  13.  Time  of  arrival  of  quenched  products  at  a  specified  density  as  a  function  of  monitoring  position  for  (a)  Case  I,  (b)  Case  II,  (c)  Case  III. 


Table  4 

Velocity  [mm/psec]  of  monitored  (specified)  density  front  observed  in  the  vacuum  region  for  each  Cases  I,  II  and  III. 


Monitoring  density 

Case  I 

Case  II 

Case  III 

1(T8  g/cm3 

19.64343 

18.99816 

19.67554 

10-7  g/cm3 

18.46885 

17.87125 

18.49775 

1CT6  g/cm3 

17.31006 

16.73676 

17.34145 

10-5  g/cm3 

15.81749 

15.27493 

15.85946 

10-4  g/cm3 

13.76884 

13.13290 

13.82020 

10-3  g/cm3 

10.50874 

8.14376 

10.50174 
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Fig.  14.  Comparison  of  the  two  different  closures  at  probe  locations:  Monitoring  TOA  of  front  of  density  10  8  g/cm3  of  (a)  Case  I  and  (b)  Case  II. 


Fig.  15.  Comparison  of  the  two  different  closures:  the  density  at  x  =  680  mm  versus  time  for  (a)  Case  I  and  (b)  Case  II. 


Table  5 

Comparison  of  the  velocity  determined  by  a  monitoring  density  for  closures, 
PT-EQB  and  PEQB  [mm/ps]  for  Case  I. 


Monitoring  density 

PT-Equilibrium 

P-Equilibrium 

10-8  g/cm3 

19.64343 

13.47709 

10-7  g/cm3 

18.46885 

12.90323 

10-6  g/cm3 

17.31006 

12.28501 

10-5  g/cm3 

15.81749 

11.54734 

10-4  g/cm3 

13.76884 

10.57082 

10-3  g/cm3 

10.50874 

8.952551 

as  a  specified  volume  ratio  to  define  the  mixture  equation  of  state.  One  must  ensure  that  the  Liu-Smoller  conditions  hold  for 
the  mixture  EOS  in  the  limiting  case  of  a  vacuum.  The  methods  illustrated  here  are  fairly  robust  and  are  easy  to  implement 
in  any  cell-based  finite  difference  code  framework.  We  feel  that  this  is  a  welcome  addition  to  a  set  of  basic  methods  in 
an  existing  simulation  code  framework  for  interactions  between  multiple  materials  that  are  subjected  to  extreme  pressure 
gradients.  In  cases  where  there  is  a  violent  expansion  of  material  into  vacuum  or  near  vacuum  regions  it  is  essential  to 
compute  sensible  and  physically  defined  fluxes,  as  defined  by  the  underlying  models  in  order  to  avoid  severe  computational 
difficulties.  Thus  we  believe  the  methodology  presented  here  will  have  widespread  applications  and  extensions  to  a  much 
larger  class  of  problems. 
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